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ABSTRACT 


Combined Numerical/ Analytical Perturbation Solutions of the Navier-Stokes Equations 
for Aerodynamic Ejector/Mixer Nozzle Flows. 

Lawrence Justin De Chant 

Texas A&M University, College Station, Texas 77843 
Principle investigators: Dr. Malcolm J. Andrews 

Dr. Jerald A. Caton 

Submitted to the Office of University Programs of 
the NASA Lewis Research Center 
in partial fulfillment of the requirements of the 

NASA Graduate Researchers Program and Grant NGT-5 1244 

In spite of rapid advances in both scalar and parallel computational tools, the large number of variables 
involved in both design and inverse problems make the use of sophisticated fluid flow models impractical. 
With this restriction, it is concluded that an important family of methods for mathematical/computational 
development are reduced or approximate fluid flow models. 

In this study a combined perturbation/numerical modeling methodology is developed which provides a 
rigorously derived family of solutions. The mathematical model is computationally more efficient than classical 
boundary layer but provides important two-dimensional information not available using quasi- 1-d approaches. 
An additional strength of the current methodology is its ability to locally predict static pressure fields in a 
manner analogous to more sophisticated parabolized Navier Stokes (PNS) formulations. To resolve singular 
behavior, the model utilizes classical analytical solution techniques. Hence, analytical methods have been 
combined with efficient numerical methods to yield an efficient hybrid fluid flow model. 

In particular, the main objective of this research has been to develop a system of analytical and numerical 
ejector/mixer nozzle models, which require minimal empirical input. A computer code, DREA Differential 
Reduced Ejector/mixer Analysis has been developed with the ability to run sufficiently fast so that it may be 
used either as a subroutine or called by an design optimization routine. Models are of direct use to the High 
Speed Civil Transport Program (a joint govemment/industry project seeking to develop an economically viable 
U.S. commercial supersonic transport vehicle) and are currently being adopted by both NASA and industry. 

Experimental validation of these models is provided by comparison to results obtained from open 
literature and Limited Exclusive Right Distribution (LERD) sources, as well as dedicated experiments 
performed at Texas A&M. These experiments have been performed using a hydraulic/gas flow analog. Results 
of comparisons of DREA computations with experimental data, which include entrainment, thrust, and local 
profile information, are overall good. Computational time studies indicate that DREA provides considerably 
more information at a lower computational cost than contemporary ejector nozzle design models. Finally, 
physical limitations of the method, deviations from experimental data, potential improvements and alternative 
formulations are described. 

This report represents closure to the NASA Graduate Researchers Program and Grant NGT-5 1244. 
Versions of the DREA code and a user’s guide may be obtained from the NASA Lewis Research Center 
grant monitor: 
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Shari Nadell 
Mail Location 77-2 
NASA Lewis Research Center 
21000 Brookpark Road 
Cleveland, OH 44135 

nadell@lerc.nasa.gov 
(216) 977-7035 

Additionally, a related, but more theoretically detailed discussion of this research is contained in the 
dissertation by the author entitled: 

Combined Numerical/Analytical Perturbation Solutions of the Navier-Stokes Equations 
for Aerodynamic Ejector/Mixer Nozzle Flows. 

Texas A&M University, May 1997 


which may be obtained: 


UMI Dissertation Services 
A Bell & Howell Company 
300 N. Zeeb Road, Ann Arbor, Michigan 48106 
1-800-521-0600 
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1. INTRODUCTION 


1.1 Objectives and summary 

In spite of the rapid advances in both scalar and parallel computational tools, the large number and breadth 
of variables involved in both design and inverse problems make the use of sophisticated and even relatively 
simple (parabolized or boundary layer) fluid flow models impractical. With this restriction, it may be 
concluded that an important family of methods for mathematical/computational development are reduced or 
approximate models. In this study a combined perturbation/numerical modeling methodology is developed 
which will provide a rigorously derived hierarchy of solutions. These solutions are characterized by varying 
levels of complexity versus analytical fidelity. Additionally, the solutions to these models utilize analytical 
solutions to resolve singular behavior as required. Hence, classical methods are to be combined with efficient 
numerical methods to yield an efficient and original class of fluid flow models. 

In particular, the main objective of this research is to develop a system of analytical and numerical 
ejector/mixer nozzle models, which require minimal empirical input. A code, DREA Differential Reduced 
Ejector/mixer Analysis has been written with the ability to run sufficiently fast such that it may be used either 
as a subroutine or called by an design optimization routine. Models are of direct use to the High Speed Civil 
Transport Program and are in the process of being adopted by both NASA and industry. Experimental 
validation of these models is provided by comparison to results obtained from the literature, industry 
proprietary data, as well as, a dedicated experiment performed at Texas A&M. These experiments have been 
performed using a hydraulic/gas flow analog model and provide information about the inviscid mixing 
interface structure. 

Since this reduced/approximate approach is applicable to a wide range of fluid flow problems, an 
additional family of test problems described in a preliminary manner, come from geophysical fluid dynamics. 
These problems directly influenced the linearization and solution methods applied to the main aerodynamic 
problem. 


1.2 Motivation 


1.2.1 Ejector nozzles 

As indicated in the summary, the basic physical problem of interest in this project is the mixing and 
operation of ejector-mixer nozzles. Although the physics of ejector nozzles is described in greater detail 
subsequently, it is useful to define the type of devices that are of interest. An ejector is a relatively simple, 
passive mixing/pumping device which serves to entrain (or pump) fluid from a secondary stream, mix it with a 
primary, high energy stream; thus obtaining a mixed (and potentially uniform) exit stream of greater mass flow. 
This process is shown schematically in Figure 1.1. 
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High speed primary jet (engine core); 

entrains fluid (viscous and local pressure differential), 

thus causing secondary stream 


FIGURE 1.1 Schematic of ejector nozzle operation. 


The concept underlying the use ejector/mixer nozzles for aerodynamic applications is driven by two 
conflicting design requirements. One is to reduce takeoff noise to an acceptable level, while the other is to 
maintain a predetermined gross thrust. The conflicting requirements of this problem are revealed by examining 
the relationship for acoustic power (a measure of noise) (Lighthill (1952) (1961), (1963) and Curie (1961)). 

P acoustic V ( 1 . 1 ) 

where N is an exponent varying between 3-8 and V is an ideal exhaust velocity. The smaller value of N 
corresponds to a fully supersonic jet, while the larger value is that of a fully subsonic jet. Though acoustic 
measurements are notoriously complex, industry and NASA designs tend to based upon die larger values of N. 
An alternative to this relationship which indicates this design relationship compares exhaust velocity versus 
sideline noise is presented in Figure 1.2. 

Also presented in Figure 1.2, is the maximum sideline noise threshold labeled FAR36 Stage 3, which is a 
regulation for the maximum ‘‘airport” noise to be accepted. Clearly, noise is strongly dependent on jet exit 
velocity. 

Thrust on the other hand may be estimated by the ideal thrust momentum relationship: 


F=mV 


( 1 . 2 ) 



HSCT Sideline Noise 

Sideline Noise @ Mach 0.32, 689 ft. 



L • • * 

1000 2000 3000 


Ideal Exhaust Velocity, ft/s 


FIGURE 1.2 Side line noise versus ideal exhaust velocity and design threshold for supersonic transport 

concept. 


Now, if thrust is to be maximized and we must reduce exit jet velocity to a minimum, the only term 
available to control is the mass flow rate of the propulsion system. Hence the objective here is to efficiently 
maximize the flow rate. Ejector nozzles, since they are an elegant passive device, have the potential for 
providing this mass flow augmentation in a very efficient manner. The desired modification of the flow field 
by use of an ejector/mixer nozzle is shown graphically in Figure 1.3. 
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HSCT Mixer/EjectorVelocityProfile 



a 


Top View 
Vel. Profile 


Velocity Profile Without Mixer/Ejector 



Velocity Profile With Mixer/Ejector 


FIGURE 1.3 General comparison of exit velocity field untreated by ejector versus treated by ejector. 


1.2.2 Installation and basic flight hardware 

From the above application, it is apparent, that the ejector nozzle is to be inserted in the engine/cycle flow path 
to entrain additional mass. As such the ejector is installed in the flow path after the engine core as shown in 
Figure 1.4. 


HSCT Propulsion System 


Noise suppression through the entrainment of 
secondary low-speed air 


vortical mixer 
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The close analogy between an engine concept fitted with an ejector and a high bypass turbofan is to be 
noted (Morrison (1996)). In fact, one might consider a turbofan to be a forced pumping, poorly mixed, i.e. 
unshrouded ejector. Also, the terminology “ejector/mixer” nozzle refers to two concepts: 

1. Ejector implies entrainment of fluid from ambient conditions. In our work the ejector nozzle refers to the 
overall system. This includes a set of inlets and chutes (shown in cross Section in Figure 1.5) and a 
“mixing duct” 

2. The mixer refers to the convoluted splitter plate system. Of course, the goal of the complex lobing is to 
enhance the mixing and entrainment of the flow 


Two-Dimensional Mixer/Ejector Design 



Vortical Mixer/Ejector Cross-section 


FIGURE 1 .5 Cross section of mixer/ejector nozzle system showing convoluted splitter plate of mixer. 


1.3 Literature Review 

Ejector/mixer nozzles have been used for industrial (mixing, jet pump) and aerodynamic (noise 
suppression, thrust augmentation) applications for many years. An example of this is an atomizer, a simple 
entrainment device used to produce a fine aerosol or mist. As such a significant body of theoretical, 
computational and experimental literature exists for these devices. Consider the bibliography of over 1600 
papers related to ejector research compiled by Porter and Squyers (1978) before 1978. A portion of this basic 
literature is described to show that intermediate level models (more powerful than control volume analyses, and 
less costly than computational simulations) are needed. Then, current models and research that have been 
developed to fill this gap are compared. 

The simplest modeling approach to ejector/mixer nozzle formulation involves using simple 1-d, or lumped 
parameter control volume formulations. An incomplete list of workers who have developed this type of 
analysis include: Presz et al. (1986) and (1990) as well as Braden et al. (1982), Alperin and Wu (1983a) and 
(1983b), Addy et al. (1981), Quinn (1973), Dutton and Carroll (1988) and (1988), Keenan et al. (1950), and 
Ginoux (1972). Alternative viewpoints, for example the lifting surface analogy of Bevilaqua (1978), have also 
been developed. Distinguished from the basic mixing analysis, which predicts ideally mixed, 1-d exit 
conditions are the supersonic ejector models (downstream flow does not influence inlet conditions), discussed 
initially by Fabri and Siestrunck (1958), which predict secondary inlet conditions for this problem. Most other 
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workers in this flow regime have followed this approach, e.g. Addy et al. (1981) and Dutton and Carroll 
(1986). Many of these workers have provided experimental measurements for basic 1-d quantities such as 
entrainment, thrust etc. Abe et al. (1992), Chandrasekhara et. al. (1991), Fernando and Menon (1993), Goebel 
and Dutton (1991), and Yang et al. (1985) provide experimental results useful for simple (planar shear layer) 
mixing flow and ejector studies. An alternative formulation, De Chant (1993), provides estimates of design 
information (geometry and operating conditions) using a similar integral approach. 

Unfortunately, none of the control volume based approaches can make any meaningful prediction about 
the streamwise length required to achieve a desired level of mixing. This length or equivalently rate of mixing 
is necessary for any aerospace application where length required for mixing translates directly to weight, a 
critical flight design constraint. To obtain this information, more complete mathematical models are needed. 
Early models employed boundary layer or 2-d, inviscid (method of characteristics) formulations to provide this 
type of information. These models require that the primary stream be supersonic. Examples of inviscid 
methods (corrected integral with boundary layer losses) include Chow and Addy (1964) and Chow and Yeh 
(1965). Boundary layer formulations include the studies by Gilbert and Hill (1973) and Hedges and Hill 
(1974a) and (1974b). These method-of-characteristics and boundary layer models provide considerably more 
information than their integral, 1-d counterparts, though again, at much greater computational cost. The work 
developed here provides a model of intermediate complexity, more powerful than the control volume 
formulations, yet simpler (and in several ways more general) than the either the inviscid computations or the 
boundary layer methods. 

The previous paragraph emphasized the need to minimize the mixing length of ejector/mixer nozzle for 
aerodynamic applications. One of the methods currently being used to improve the mixing performance of 
ejector mixer nozzles involves replacing the straight splitter plate by a convoluted or “chuted” mixer plate. 
This has the effect of inducing streamwise vorticity in the flow and accelerating the mixing rate. Here also, a 
large number of experimental and computational studies have been performed. Since the flow is rather 
complex, full Navier-Stokes simulations have been typically used to simulate these flows. Again a partial list 
of experimental and full simulations include: Presz et al. (1986), (1987a), (1987b), Skebe et al. (1988a), Lord 
et al. (1990), Tillman et al. (1992), Povinelli et al. (1980), Malecki and Lord (1990), Debonis (1992), Barber 
and Anderson (1991), Barber et al. (1988), Keith et al. (1993), Abdol-Hamid et al. (1993), Booher et al. 
(1993), Qui (1992), Elliott et al. (1992), Tew (1992) and O’Sullivan (1993). Strictly experimental results 
include: Paterson (1984), Werle and Paterson (1987), Skebe et al. (1988b), Koutomos and McGuirk (1989), 
Manning (1991) and McCormick (1992). Although the simulations and, by their very nature, experiments are 
potentially excellent representations of full scale ejector mixers, the results are far too expensive in terms of 
cost, manpower and expertise to interact directly with design or inverse design computations. 

Alternatives to the previous research methods, include this study, a code developed at MIT, Fung (1995), 
as well a code developed by Boeing and Clark (1995) are compared in Table 1.1. As can be seen from Table 
1.1, the DREA code provides a significant improvement over the other current models, all of which have been 
developed to fill the need for a fast, yet sufficiently powerful computational model. 
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Criterion 

DREA (Texas A&M) 

Fung (MIT) 


Vortical Mixing 

yes; analytically based 
turbulence model 

yes, piecewise model and 
scaling functions 

no, multiple lobe 
only 

governing equations 



l-d 

primitive variable 
profiles 

yes 

no 

two stream 
model 

empiricism 

limited 

moderate (“scaling” 
constants) 


compressible 

expressly developed for 
compressible flow 

1 -d and extensions 

1 -d and 
extensions 

internal shock /choke 
structure 

yes; Fabri-choke and 
back pressure limited 

no, pumping only 

current version 
no 

computational speed 
(suitable for 
preliminary 
design/optimization) 

fast (user may choose 
control volume and/or 
mixing) 

fast 

fast 

experimental validation 

yes, literature, 
proprietary data and 
experiments 

yes, literature and 
experiments 

yes, experiments 

publications 

journal articles, 
conference papers, 
dissertation and reports 

thesis, conference papers 

conference 
papers 
(main code 
proprietary) 

status availability 

preliminary summer 
1996, fully tested and 
documented 05-97 

complete 

complete 


TABLE 1.1 Comparison of ejector-mixer models suitable for preliminary design/optimization. 


Analytical solution methods, numerical differencing and solution algorithms, and formal perturbation 
methods will be used thorough out this project. No attempt is made to give a formal historical literature 
development of these methods since they are truly subjects unto themselves. However, the fundamental and 
relevant references are provided as we develop our specific problem. 

1.4 Mathematical framework 

This section describes the overall mathematical framework for the ejector problem. The requirements for 
the analysis are presented first. Then the basic perturbation hierarchy, which is used to develop the systems of 
governing equations, is described. 

1.4.1 Overall analysis 

The model developed here describes an analysis code, i.e. one for which geometry and operating 
conditions have been specified. Performance characteristics (thrust, noise, degree of mixing, and pumping) are 
predicted in a very efficient manner by this code. Ultimately, this type of tool may be used to generate off line 
performance information in the form of ’’ejector" maps, which are tabular list of performance characteristics. It 
is also sufficiently computationally efficient to be run in a "real time" fashion. Table 1.2 presents the 
requirements of this analysis follow: 
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Computational 

efficiency 

Hierarchy of 
complexity 

Minimum of 
empiricism 

Experimentally 

validated 

Documentation 

Reduced equation 
formulation 

Formal 

perturbation 

First principle 
methods 

Literature and 

dedicated 

experiment 

journal, dissertation, 
reports 


TABLE 1 .2 Requirements for computational simulation tool. 

The efficiency of this method also permits the realistic use of this analysis code as the "base solver" for an 
iterative, inverse design analysis. In this problem, geometric and operating conditions are parametrically 
adjusted to meet a series of overall design constraints. This inverse design methodology can be quite efficient, 
especially if a reasonable first iteration estimate is available. A simplified method that may be used to estimate 
operating and geometric parameters has been developed by De Chant (1993). Hence, there is an overall family 
of codes that provide a useful preliminary design capability for the design and implementation of ejector 
nozzles. 

1.4.2 Formal perturbation hierarchy 

As implied in the previous section, one of the main strengths of the overall modeling methodology is 
development of a formal perturbation hierarchy. This type of expansion method permits us to "apriori" 
estimate the error associated with a particular approximation, as well as, formulate corrections to the 
approximation. Of course, approximate methods in fluid mechanic have been well developed over the years. 
For example considering steady, viscous problems only, consider the hierarchy Anderson et al. (1984) and 
Anderson (1989): 


Perturbation Method/Approximation 

Comments: 

Integral/Control volume 
(.i.e. 1-d or quasi-l-d) 

Satisfy integral forms of conservation equations. Cannot 
predict local values. 

Current DREA model 
0(1) equation system 

Parabolic in terms of conservative flux quantities, primitive 
quantities: u, v, p, p, M, locally predicted, PNS subsonic 
pressure approximation 

Boundary layer 

Parabolic in terms of u, v for all flows; p(x), thin layer, p(x) 
must be specified by an auxiliary relationship 

Viscous-shock layer 

Parabolic for supersonic flow; u, v, p(x) specified by auxiliary 
relationship, p(x,y) specified by inviscid y-momentum; PNS 
(Vigneron et al (1978)) subsonic pressure approximation. 

Parabolized Navier-Stokes 

Parabolic (in space) for supersonic flow, elliptic subsonic 
regions; u, v, p(x), p(x,y) Viscous y-momentum. PNS subsonic 
pressure approximation 

Full Navier-Stokes Analysis 

Parabolic-hyperbolic, fully viscous equations. Always solved 
using unsteady or iterative method. 


TABLE 1 .3 Comparison of classical and current approximation methods in fluid dynamics. 
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Each level has certain strengths and weaknesses associated with it. Complexity and cost increase at an 
extraordinary rate as we work our way from down Table 1 .3. 

The work described in this section is to develop a system of partial differential equations describing 
viscous mixing. Expansions that are less complex than the boundary layer family but provide more information 
(local mixing field 1 ) than the integral methods are of special interest in this study. Like a quasi- 1-d model, 
though, it is desirable to retain (at least approximately) the local prediction of static pressure , rather than the 
external imposition of a pressure field, predicted either using free stream information or through a global mass 
conservation constraint, as would be done for a classical boundary layer analysis. As is shown in Section 2, 
this local pressure prediction requirement adds considerably to the complexity of the model. 

The fundamental governing equations for this project (lowest order system) are developed in this section 
and higher order expressions are deferred to Appendix A. 

To this end, the development begins by considering compressible, two-dimensional Reynolds equations 
Anderson et al. (1984) or White (1991). Note that the turbulent Prandtl number has been assumed to be close 
to unity, hence dissipation terms in the energy equation, may be neglected. Due to the high Reynolds number 
asymptotic nature of these equations, the molecular viscosity and the molecular Reynolds numbers are ignored 
(See Schlichting (1979) or Tennekes and Lumley (1972)). The resulting equations are: 

x-momentum: 

4-(P^)+4r(f m > + <U) 

dx dy dx dy ox 

y-m omen turn: 

— (puv ) + ^-(P v 1 ) + f£ - • 4 - (P“' v '> -4r(t»' v '> <’ - 4) 

dx dy dy dx dy 

energy: 

<>•’> 

and continuity: 


-^-(pu) + -^-(pv) = 0 0-6) 

dx dy 


To proceed with the analysis, it is recognized that when non-dimensionalized, there are a number of small 
terms” inherent to the governing equations. They come from both the magnitude of the dependent and 
independent variables. The flux difference between streams is relatively small, as is the mixing layer thickness 
compared to the streamwise distance. Formally stated this yields: 
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(1.7) 


s _ (U io~ Uia) K £ _ (Gw Gy) _ Smix 

(U\ 0 +U 20 ) G (G I0 + G 20 ) L 


where 8^ refers to the mixing layer thickness and G=pu 2 +p. It should be noted that e«ec may not be small for 
all cases. Indeed s«e G =1 in the free jet limit. In spite of this, the governing equations will be shown to 
provide reasonable results despite the fact they may be “formally tenuous” approximations. Support for 
effectiveness of the free jet limiting case may be found in Schlichting (1979). 

Though it has not explicitly stated at this point, this expansion parameter e which is valid for turbulent 
flow will have a form analogous to that of Ting (1959), who studied laminar mixing between two streams. Of 
course the essential stretching transformation, as well as, the Poincare’ expansion methods are fundamental to 
virtually all perturbation (especially singular perturbation) approximations. Nayfeh (1973), Kevorkian and Cole 
(1981), Van Dyke (1975) and Cole and Cook (1986) provide developments for a wide range of classical 
problems where regular and singular perturbation methods have been applied. 

It may be noted, that the finite or bounded form of internal ejector flow makes the classical definition and 
then “patching” of inner and outer solutions problematic for the problem developed here. As such, stretching 
transformations are used to develop a solution valid through out the ejector flow field. 

Since it is necessary to derive differential equations that can "see" or resolve the mixing layer itself, it will 
be necessary to rescale (stretch) the cross-stream or "y" coordinate: 

n/i\ - * 5 / d 

x * 0(1) y = siy — = — — ^ ( 1 . 8 ) 

dy £l dy 


This is, of course, the classical boundary layer scaling as introduced by Van Dyke (1975). The selection to 
rescale the cross-stream variable to permit resolution of the mixing layer is the connection between the partial 
differential equation based models and the simple control volume based models. Literally if it is assumed that 
for the integral analysis that the cross-stream scale is 0(1) all "viscous" terms on the right hand side are lost, 
governing equations revert to the control volume formulations. 

Proceeding with the scaled problem, we further estimate the magnitude of the Reynolds stress terms. This 
starts, by considering the Boussinesque approximation (introduction of an effective turbulent viscosity 
(Tennekes and Lumley (1972), Hinze (1959) or Rodi (1993)) for the streamwise momentum flux: (pu 2 )V 
(Schlichting (1979) and Tennekes and Lumley (1972)): 


GV = (puu)'v ' « ~v ej j' G 


i dy 


~.(pu 3% ) 


(1.9) 


This quantity is chosen since the streamwise momentum flux is a conservative value, rather than simply 
velocity, which is a primitive variable. This choice of variables helps to overcome one of the main objections 
raised to using mixing length based algebraic turbulence models, which do not even conserve momentum 
(Tennekes and Lumley (1972)). To apply this to our problem, the momentum flux transport must be must 
related to the transport of velocity fluctuations by turbulence, i.e. puV by introducing the definition: 
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( 1 . 10 ) 


pu V = 


(p uuyv_ _ V JJ_Q_^G d 
u u \ dy 


~.(pu ) 


This is considered to be an acceptable approximation, since the effective viscosity closure is a definition 
already. To apply this equation, we must linearize terms like: 


1 _ 

1 

u 

1 

• * o 

U 

u 0 + w, s l 


_ _ 


( 1 . 11 ) 


where U=1/2 (Ui 0 +U 2 o) Since, u 0 * is our 0(1) approximation, it will be acceptable, to 0(1), to write: 


u 


U 


1 -u.s 2 +. 


( 1 . 12 ) 


Hence our 0(1) approximation for the Reynolds stress must be: 


pu'v' 


eff 


yv ) 


Os' dy 


-.(pu 2 ')+. 


(1.13) 


The physical basis for this scaling and the turbulence model values, such as, v efl /U, are discussed in Section 
2.5. 


Using equation (1.13) and invoking the additional well founded approximation that the turbulent Reynolds 
number Pr,=0(l), the Reynolds transport of total enthalpy is written: 


pH'V 


Pr, 


v cff 

U , 


dy 


-(puH) 


(1.14) 


At this point the perturbation expansions for the various variables are introduced. The streamwise variables are 
simply expanded: 


Pf' = Pf'o + pf]e 2 + pf]e+... 

g=g0 + g,£ L 2 + S2 £+ -" g ' 


/-«*. H 


= p‘ , u , H ' 


(1.15) 
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The expansion of the cross-steam flux and velocity take more care. It is expected, that for our basically 
parallel flow that the cross stream velocity, v, is very small 0(e). The cross stream mass flux, pv, is also small, 
but to conserve mass is an order of magnitude larger 0(e ,/2 ). With these estimates, the cross-stream mass flux 
and velocity are written: 


PV - p Vo £2 + Pv) £+... 

(1.16) 

3 

V*-V o£+ V*£ 2 +V2£ 2+ — 

Using these expansions and estimates available, it is possible to scale the dependent variables, substitute 
the expansions and collect terms. Performing these operations, the two lowest order systems are obtained. The 
lowest order, 0(1), system (which is of most particular interest) is a simple set of parabolic "mixing” equations. 
Most of the analysis in the dissertation will focus on the analysis of these equations. The second system will 
be recognized to be a linear formulation of the classical boundary layer equations. See Appendix A. Hence as 
desired, this expansion has yielded a set of governing equations of lessor complexity than even the classical 
boundary equations: 

x-momentum: 


— (Puouo + p'o) 
ox 




G 


£ dy * 


(PUoUq)] 


y-momentum: 


energy: 



= 0 


— (PuoH'o) 
ox 


r[ 


1 


eff 


dy Pr, \ U J e dy 


7 (puo Ho)] 


(1.17) 


(1.18) 


(1.19) 


and mass conservation: 


d . d 

—(pu o ) + —7(pv o ) = 0 
ox dy 


( 1 . 20 ) 


This is the 0(1) system, that is solved. As it is written, it is not possible proceed, since the last term in 
equation (1.20) is not closed. To solve this system it is necessary to approximate the cross-steam mass flux 
term: 
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( 1 . 21 ) 


pvo£ 


2 ~ 


Sc, 


¥ jff 

U 


°s* Sy 


~,(pUo) 


Which is merely a linearization of our Reynolds stress closure. 

One more convenient change may be made. Since the cross-stream pressure gradient 9p/5y«0(e) l/2 , we 
can add this term to the right hand side of the streamwise momentum equation thus yielding: 


~(Puo ul + p'o) 
ox dy 


jL 

Vu, 


£ dy 


-;( Puouo + P o)] 


( 1 . 22 ) 


Finally, a state equation is included, which for our uses is an ideal (and thermally perfect) gas law. Combining 
this with the definition of total enthalpy and retaining first order terms yields: 


Po = 


r-i 

r 


pH, --pu 0 u o 


(1.23) 


Though, the solution is actually formed in terms of alternative variables it is desirable to perform an equation 
count shown in Table 1.4. 
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Equations 

Unknowns 

momentum, equation (1.17 or 1.22) 

Uo 

energy; equation (1.19) 

Ho or T 0 

mass; equation (1.20) 

Po 

mixing relationship closure; equation (1.21) 

Vo 

state; equation (1.23) 

Po 

Total=5 

Total=5 


TABLE 1.4 Equation versus unknown count for 0(1) equation system. 


It is worth noting, that this solution set, equations (1.17)-(1.22), is not merely linearized boundary layer, 
since the static pressure p 0 participates directly in the solution procedure. In a classical internal boundary layer 
solution, on the other hand, the pressure is strictly a function of x, i.e. p=p(x). This permits the estimation of 
p(x) by demanding that global mass conservation be observed, i.e. total mass= jpudy. Our method, 
alternatively, is more closely related to local, 1-d and quasi- 1-d compressible modeling relationships. The 
benefit to us, is the maintenance of a strong thermodynamic coupling between variables and the ability to 
predict p(x,y) on a local basis. That on a local basis p(x,y) may show strong variations both streamwise and 
across the mixing layer, seems to contradict the y-momentum equation. To satisfy this equation it is formally 
required to model equations for matching inlet pressure. In reality, though, we will often need to model 
problems with large static pressure imbalances. Apparently, due to the local pressure computations, our 
equations still provide a useful model even in these situations. Classical boundary layer could never achieve 
this, since p=p(x) only. Other quantities, such as, Mach number and total pressure are supplemented by their 
own specific definitions. 

It has been stated that equations (1.17)-(1.22) are not a straight forward linearization of the boundary layer 
equations. The basis for this statement comes from the fact, that linearizations were developed using an 
alternative transformation (the use of the conservative fluxes, pu 2 +p, puH and pu) rather than primitive 
variables. Other examples of the use of a transformation and subsequent linearization for incompressible and 
compressible boundary layer flows is the classical von Mises Transformation (Schlicting (1979), Von Karman 
and Tsein, (1938) Schetz and Jannone (1965)). We remark that Schetz and Jannone (1965) compare primitive 
variable linearization versus von Mises transformation and linearization and conclude the latter method is 
superior for a wide range of flows. This studiy helped to motivate our use of the conservative flux 
transformation and linearization employed in this study. 

Before closing the perturbation development is worth describing the system of equations consistent with 
treating the cross-stream or ”y" variable as 0(1). As indicated earlier, they take the form: 
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d 

dx 


( pu 0 Uo + Po) = 0 


—(p 0 u 0 Hl) = 0 
ox 


(1.24) 


d_ 

dx 


(pul) = 0 


which is precisely what is to be expected. Integration of these relationships provides no information as to the 
streamwise length scale. These are literally control volume formulations. Though control volume formulations 
provide no streamwise mixing length scale, we will use control volume and integral methods to estimate initial 
conditions for the mixing portion of the code. These analyses are described in detail in section. This completes 
the formal perturbation development of the systems of equations used in this project 


1.5 Basic physical problem 

This section describes the basic physical arrangement of the ejector/mixer problem. Though ejectors have 
already been described briefly, it is the goal of this section to concentrate upon the physics and variety of flow 
characteristics. An ejector is a relatively simple, passive mixing/pumping device that serves to entrain (or 
pump) fluid from a secondary stream, mix it with a primary, high energy stream; thus obtaining a mixed (and 
potentially uniform) exit stream of greater mass flow. The additional mass entrainment is cause by two major 
effects: (1) inviscid pressure imbalance, which literally sucks fluid into the stream from the secondary inlet and 
(2) an effective viscous component, which drags fluid in from the secondary inlet. 

It is necessary to be somewhat careful in the use of these two end member decompositions, in that through 
the momentum equation, pressure and effective viscous terms are coupled. Therefore, in reality, no simple 
decomposition is appropriate. This coupling of mixing and pressure will become apparent as the pumping 
computation methodology is developed in detail. Additionally, we need to define ‘effective viscous 
components of the entrainment as terms which are caused by turbulent mixing effects. These effects are 
typically estimated using the concept of a turbulent viscosity, but are really strong, non-linear fluctuation type 
terms, i.e. u’v’ terms and have no relationship to molecular viscosity. Though molecular viscosity is certainly 
present in any physical flow, the Reynolds number is assumed to be large enough such that molecular viscosity 
effects are wholly negligible, since they are multiplied by terms like l/Re mo |., visc . 

An interesting limiting case of the pumping computation is the free 2-d jet. This would be the case for an 
ejector where the area ratio A 2 /A, -> oo. In this limit, the pressure balancing control volume concept that 
would be applied is not suitable. There the entrainment is truly due to local effective viscous (turbulent) 
interaction. Following Schlichting (1979) the local self similar velocity relationship for a two-dimensional 
turbulent jet is quoted: 


V3fHVl t on u2 r,-sr y 

u = — — 1 — tann rj) n = a— 

2 V x ) v ' x 


(1.25) 


where k is the net momentum flux, ju 2 dy and c is an empirical turbulence constant, i.e. 1/a ac v eff . Now, to 
compute the mass flow rate entrained we integrate over the flow field, to obtain: 


35 



( 1 . 26 ) 


Q = 


_ V3 f kx 


2 V <7 


1 

2 “ 


J(l - tanh 2 



From this equation, it is apparent that entrainment increases as we proceed from the origin. To our 
dissatisfaction, it is noted that this flow rate integral has entirely neglected the contribution of the primary 
stream! This failure is due to the fact, that for our case with its infinite secondary area, the flow rate from the 
primary is negligible. This emphasizes the limitation of the free jet approximation for use in finite ejectors. 
Ejectors are finite devices. Integral quantities may not be estimated using infinite space methods. Hence the 
proper use control volume methods is to predict entrainment and thrust. As a side note, though, it may be 
noted that locally non-integral quantities such as velocity etc. can be acceptably modeled using infinite method 
if we are far from confining walls. 

Depending upon the operation of the ejector net thrust may also increase. As indicated previously, it is the 
ability to maintain (or even improve) thrust, while reducing the exit velocity which is fundamental to the 
application for high speed civil transport aircraft. 

Since ejectors involve the mixing of two streams of fluid which (for compressible gas flows) may be either 
supersonic or subsonic, several possible flow types or flow regimes may exist. These flow regimes are strongly 
characterized by the extent of supersonic or subsonic flow. This is directly related to the physics of subsonic 
versus supersonic flow. 

Supersonic flow is properly modeled using a parabolic-hyperbolic equation set. The need for a hyperbolic 
system stems from the fact, that due to the supersonic nature of the flow, information or disturbances, are 
convected downstream more rapidly than can be transmitted upstream by molecular interaction (i.e. speed of 
sound). Hence, there is no chance of sending any information upstream. As one would expect form the 
mathematical classification, (parabolic-hyperbolic) flows of this type are dependent solely on their initial 
conditions. In contrast, the convective speed of a subsonic flow field is less that the molecular signal 
propagation velocity, i.e. speed of sound. As such, downstream signals can eventually propagate back 
upstream. Here, the flow field and even the initial conditions depend upon the downstream conditions. The 
modeling equations for this case are parabolic-elliptic. 

Ultimately simplified equations that have a single character, i.e. parabolic are developed. We will, 
however, supplement this parabolic system, by initial conditions which (depending on the flow regime) will 
respect the potential for downstream dependence or downstream independence. It is apparent, that all of the 
problems of interest will virtually always contain some region of subsonic flow; forcing us to respect some 
form of a downstream constraint. 

The first flow that is considered is a supersonic inlet, which has a sufficiently large back pressure (i.e. 
ambient pressure) to cause the supersonic primary stream to go through a family of oblique shocks, ultimately 
terminating in a strong normal shock. This is illustrated in Figure 1.6. In Section 2 an approximate criterion for 
the critical exit pressure needed to cause this shock is developed. With this “shock down” the exit flow is fully 
subsonic. As such, the local pressure must match the external giving a constraint that is used to estimate the 
secondary entrainment. This is also developed in detail in Section 2. 
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M <1 

secondary 



Primary stream shock down to subsonic flow 
terminating in normal shock 


FIGURE 1 .6 Ejector nozzle in subsonic (back pressure dependent mode) with shock down of primary 

stream. 


This shock down criterion is dependent upon complete mixing between the two streams. The complete 
mixing assumption is reasonable for a long ejector, but unrealistic for many short shroud ejectors. Often, 
due to insufficient mixing, two distinct streams exit the ejector: one supersonic and one subsonic. See 
Figure 1.7. For this type of flow, though, the secondary entrainment is dominated by the exit pressure 
because it is subsonic. When this occurs, the entrainment should be predicted, as before, on the basis of the 
exit pressure. 


M <1 

secondary 



M . >1 

primary 


M <1 

2,exit 


M >1 

l,exit 


FIGURE 1.7 Ejector nozzle in back pressure dependent mode due to poorly mixed secondary stream. 


Finally, for sufficiently low back pressure, the flow is fully supersonic at the exit plane. This situation may 
occur if the primary stream accelerates (expands) while the secondary stream (also expands) but chokes. As 
shown in Figure 1 .8 this expansion/choking phenomenon causes the streamlines separating the two flows to 
form an aerodynamic or Fabri (Fabri and Seistrunk (1958)) choke. Clearly then the exit stream is supersonic, 
and as such, is independent of the back pressure. The local effect of the subsonic stream does, though, 
influence the secondary entrainment. The information, though, that is sent into the secondary inlet is no longer 
a pressure constraint, but the fact that the flow has choked. An entrainment model based upon this scenario is 
developed in Section 2. 
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M=1 

Aerodynamic, "Fabri" choke 



FIGURE 1.8 Ejector nozzle in back pressure independent mode, exhibiting aerodynamic or Fabri choke. 


Finally a special case of the aerodynamic choking phenomenon occurs when the choke forms within the 
secondary inlet itself. This situation is called saturated supersonic flow and is represented by Figure 1.9. 


M=T 

Aerodynamic, "Fabri" choke 
(special case of supersonic ejector flow) 



M >1 

primary 


M >1 

exit 

P<P 

crit 


FIGURE 1 .9 Ejector nozzle in back pressure independent mode. Aerochoke in secondary inlet causing 

classical “supersonic saturated” operation. 


As one would expect, this flow is completely specified by upstream conditions, with no chance of any 
downstream influence. 

Due to the possibility of mixed or transonic flow fields, ejector nozzles have a rich physical and 
mathematical character. In Section 2 analytical and numerical models to deal with this varied character are 
developed. In Section 5 long-term solutions to problems, which may be directly traced to the multiple 
regime flow, fields that we have just discussed are developed. 
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2. THEORETICAL DEVELOPMENT 


2.1 Overall solution methodology 

The modeling algorithm developed in this study is conveniently divided into several modules or 
components. These modules are combined to approximate the flow structure of the ejector. Although the 
ejector is an integrated device (secondary stream inlet static conditions are flow dependent for any degree of 
subsonic flow, i.e. elliptic flow, in the secondary inlet stream) it is far more efficient to model this dependency 
in a decoupled fashion. This decomposition is acceptable, since the secondary stream conditions are in a large 
sense governed by inviscid phenomena. Hence, the secondary stream pumping may be estimated using ideal 
relationships, (though we have introduced modifications to account for mixing losses and the effects of 
incomplete mixing on the initial conditions). These methods are described in Section 2.2. The actual turbulent 
mixing region is model using a single pass, parabolic marching method, which provides flow field/mixing 
information and a more complete estimate of two and three dimensional effects within the mixing section of the 
ejector and are developed in Section 2.4. Consider the "flow chart" description of these modules and this 
iterative solution. Additionally, the possibility of an inverse design iterative loop is shown, though this level of 
iteration is not within the scope of the current project. 
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♦Part of the data reported in this section is reprinted with permission from Interface Wavelength between supersonic 
jets and subsonic flowfields. by De Chant, L. J. Seidel, J. A. and Andrews, M. J. 1996. AIAA Journal. 34, 1946-1948. 
Copyright 1996 by American Institute of Aeronautics and Astronautics. 




FIGURE 2. 1 Flow chart of overall computational method. 
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As indicated by the flow chart. Figure 2.1, estimates of inlet conditions are principally, for the 
secondary stream flow. Further, since the parabolic marching technique which models the mixing portion 
of the flow is steady (single pass marching), but typically must model transonic flow fields, the extent of 
regions of supersonic or subsonic flow must be estimated. This task is performed using modifications of 
classical, inviscid slipline analyses, and is described in the Section 2.3. Then the mathematics and numerics 
of the rest of the flow chart are developed, which is principally the mixing algorithm portion and the 
associated turbulence model, in the remainder of Section 2 and Section 3. 


2.2 1-d Control volume analyses 

In this section a classical, 1-d, compressible, control volume based, mixer/ejector analysis is 
implemented. The basic mixing analysis consists of 1-d, conservation equations with upstream conditions 
fully specified and "mixed" or downstream conditions to be computed. Although simple, this model 
provides useful insight into the mixing process. However, the ejector mode of operation, in which the 
secondary stream static conditions are specified by downstream conditions, is not as clearly defined. Two 
possible modeling techniques for the ejector mode of operation are discussed and compared: Fabri choke 
analysis, (Fabri and Siestrunck (1958) or Addy et al. (1981)) which are back pressure independent flows 
and the use of a subsonic downstream pressure condition for back pressure dependent flow (Addy et al. 
(1981)). As implied by their descriptions, the choice of the proper modeling technique is determined by the 
back pressure conditions. 


2.2 A Introduction to control volume methodologies 

One-dimensional or control volume mixing analyses (Alperin and Wu, (1983), Addy et al. (1981), 
Dutton and Carroll, (1986) and (1988)) provide simple, but relatively accurate predictions of mixer/ejector 
operation. However, the basic integral formulations do not permit prediction of non-ideal behavior, nor do 
they provide any information concerning the rate of streamwise mixing or streamwise length scales. This, of 
course, provides the impetus for the mixing flow modeling portion of this code. Although the integral 
relationships are contained implicitly in the differential method, it is convenient to have access to these 
solutions. In fact, the control volume solutions may be seen as a "zeroth" order solution or boundary 
conditions for the higher order models and, more importantly, provide a method to predict the initial 
conditions for the partial differential equations based mixing model. 

This section seeks to describe the basics of 1-d mixing analyses. The ejector mode of operation, and its 
associated closure methods, are then discussed. The analyses associated with these closure methods, or 
boundary conditions, are derived. Finally, the analysis is verified by comparisons to available experimental 
data. 


2.2.2 Mixer/ejector formulation 

Since ejector flow modeling may be considered an extension of simple mixing, we begin by considering 1-d 
control volume conservation equations. See Figure 2.2. A number of assumptions are implicit in these 
relationships: 

(1) The entrance and exit flow fields are locally one-dimensional. When describing the exit location, this 
is a "fully mixed" assumption. Although, this is the most restrictive assumption of this analysis, it is 
shown that this model still provides good comparison to experimental data. This assumption is relaxed 
in Section 2.2.6. 

(2) Thermally perfect gases are assumed throughout, as well as, ideal gas behavior. 
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With the application of these restrictions, the conservation equations may be written: 


my + m 2 = ril) = w 


( 2 . 1 ) 


rhiui + p,Ai + m 2 U 2 + P 2 A 2 + ft^pdA = m 3 u 3 + p 3 A 3 = P mom (2.2) 


and 


miToi + m 2 To 2 ~ m 3 To 3 = E 


(2.3) 


where the subscripts, " 1 ", "2" and "3" refer to the primary, secondary and exit locations, respectively. 



The integral term in the momentum equation (2.2), is closed in an approximate manner by the relationship: 


J 


A 3 

Aj+ A + 2 


pdA » 


1 

2 [P 3 + (P, a i + P 2 M) / (A, + A 2 )][ A,+ A 2 -A 3 J 


(2.4) 


The implications and restrictions of this approximate closure are discussed in, De Chant (1993). 

For simple mixing, all upstream quantities are specified. Thus with addition of state, and constant heat 
capacity thermodynamic definitions the above equations may be solved for the downstream unknowns. In 
particular, the above relationships in terms of the mixed Mach number, for M 3 : 

[(1 -G / 2)/ + y^] M/ + 1)-Gy] M/ + ^- ( Al+A2 + if = 0 (2.5) 

2 A 3 4 A 3 


where 
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2 


( 2 . 6 ) 



REw 


Equation (2.5), has two roots, in terms of M 3 2 . One root corresponds to a supersonic solution, while the 
other root denotes subsonic flow that would result if the supersonic flow passed through a normal shock. 
We note that with specification of the exit Mach number, M 3 , all other static and stagnation quantities are 
available to us by back substitution. 

The solution in equation (2.5), also provides the essential formulation for ejector flows, in which all 
conditions for the primary stream are specified, while only the total conditions, P 02 and T 02 , are specified for 
the secondary stream. To apply the simple mixing analysis all the conditions in the secondary stream must be 
specified. This is accomplished by invoking a downstream constraint that implicitly defines the secondary 
stream static quantities. A discussion concerning the physical and mathematical modeling of this constraint is 
provided here. 

The modeling of ejector flows may be made more concrete by counting equations and unknowns as 
presented in Table 2.1. 


Unknowns 

equations 

secondary Mach number, M 2 

momentum conservation 

exit Mach number, M 3 

energy conservation 

exit velocity, u 3 

mass conservation 

exit pressure, p 3 

exit constraint 

exit density, p 3 

Mach number definition 

exit temperature, T 3 

state 

total=6 

total=6 


TABLE 2.1 Control volume analysis equation versus unknown counts. 


The specification of the downstream condition or constraint requires discussion of the physics of ejector flow. 
Three possible cases or flow regimes may be identified. These cases are (i) "mixed" flow or subsonic flow, (ii) 
supersonic flow and, (iii) a subset of supersonic flow denoted supersonic saturated flow. These regimes are 
considered in detail: 

(i) "mixed" or subsonic flow. In this case, the secondaiy stream is fully subsonic. The appropriate 
downstream boundary condition or constraint is: 

Ps = P„« s P. 

This is a classical condition for 1-d, subsonic nozzle flow. The actual computations are performed by 
guessing the secondaiy Mach number, computing the downstream pressure, and correcting M 2 until 
the downstream pressure satisfies the above constraint. 


(ii) supersonic flow. As the name implies, at the fully mixed flow, exit location, the flow is 
supersonic. Thus, the secondaiy stream must accelerate from subsonic flow to supersonic flow. The 
transonic point in this flow is often termed an "aerodynamic throat". Closure is obtained by noting 
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that the flow upstream of the aerodynamic choke is dominated by inviscid phenomenon. This 
provides the basis for the Fabri choke formulation, which is discussed in detail later. We note that 
this constraint leads to under-expanded or over-expanded nozzle operation. 

(iii) saturated supersonic flow. This regime denotes a choked secondary stream, with the choke 
occurring in the "mouth” of the secondary stream, at location (2). For this flow, specification of the 
secondary Mach number by a downstream constraint is trivial, since (M 2 =l). 

The above discussion shows that supersonic flow in an ejector requires further consideration. Before 
proceeding, it is worthwhile to note that qualitatively, subsonic (mixed) flow regimes may be expected for high 
static back pressure and supersonic or even saturated supersonic flow for lower back pressures. In practical 
ejector operation, flight from sea level to high altitude with ejector deployed is possible, as such, it is necessary 
to be able to model both back pressure dependent and back pressure independent flow regimes. 


2.2.3 Back pressure dominated flow 

As indicated previously back pressure dependent flow is characterized by a fully subsonic secondary 
stream. This is in accord with our understanding of steady supersonic/subsonic flow, since only for a subsonic 
flow is it possible for a downstream signal (back pressure level) to be sent upstream to modify the secondary 
inlet flow static conditions. Within the assumptions of classical one-dimensional control volume modeling, it is 
clear that the ejector closure for this problem is the downstream boundary condition: Ps^Pcxit^amb. 

2.2.4 Back pressure independent flow (Fabri-choke) 

Two possible supersonic flow regimes have been discussed supersonic flow and supersonic saturated flow. 
Since supersonic saturated flow is a relatively simple mixing process, it is not discussed further. Instead, the 
closure for supersonic flow, the Fabri or aerodynamic choke analysis is presented. 

Consider Figure 2.3 illustrating the development and location of an aerodynamic choke. 



FIGURE 2.3 Aerodynamic/Fabri choke analysis definitions. 


Aerodynamic or Fabri choking analysis has been discussed by several workers, Fabri, (1958) and Addy et 
al. (1981). The analysis involves invoking a new control surface located at the secondary stream aerodynamic 
choking point. See Figure 2.3. Additional assumptions include: the streams remain separate, so there is no 
mixing and the flow is isentropic for each stream. Neither assumption is exactly true, but assuming the choke 
occurs relatively close to the primary and secondary inlets and that the flow is primarily inviscid, this is an 
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appropriate analysis. These assumptions may be poor for a forced ejector where mixing enhancement is 
desired. To overcome some of these limitations, we have developed a loss mechanism (in terms of total 
pressure and total temperature) to account for the interaction between the two streams. This loss model closure 
is also presented in Section 2.2.6. 

Considering the modeling in detail and with reference to Figure 2.3, four (4) unknown quantities may be 
identified: 

A2s, A/s, M 2 , M/ s 


The related equations are: 


A 

A 2s 


M : 


L y + i 


( 1 + 


hi 

2 


M 2 2 ) 


2(7-‘) 


(2.7) 


A 

A is 


M u 

M, 


Y-l 2 

l + ^—Ml 


y-l 2 
l + h-M/l 


r+1 

2(r-0 


( 2 . 8 ) 


the so called Mach number area relationships (isentropic mass conservation), (Anderson (1985)). The equation 
set is completed with a constant area ejector cross-sectional area definition: 


A 2s + A h = A} = const. 


(2.9) 


and the isentropic momentum equation. 


PiAi[i+y Ml] + p 2 A2[i + r m/J- 


P 2s A2s[1 + Y M 2 I ] + Pis Ais[l + r Mis 2 ] 


( 2 . 10 ) 


The static pressure may be related to the total pressure using: 

y-l , ' 

P = P 0 [l + L yM 2 ]7 r r 


( 2 . 11 ) 


These (4) equations may be solved to provide a direct value for the secondary stream inlet Mach number M 2 . 
Thus, with the secondary static quantity, M 2 , specified the flow now acts like other mixing flows, where all 
quantities upstream are specified and the downstream (mixed) values are computed. 
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Equation (2.9) has implied that the ejector is locally constant area, with this area based upon the inlet 
values. This statement is equivalent to restricting the location of the aerodynamic choke close to the inlet 
plane. This is a reasonable first approximation for two reasons: 

(1) The cross-sectional area for typical designs is approximately constant, i.e. Mixer Area Ratio, 
MAR^l. 

(2) The Fabri analysis as developed here is locally isentropic. This would not be true for a long 
mixing section. Ultimately methods for estimating the aerodynamic choke location in the next section, 
as well as the losses incurred reaching that point are discussed. 

Clearly, to effectively validate this type of analysis comparison to experimental data is necessary. To start, 
we begin with both constant area and variable area, supersonic ejector experiments from Addy et al. (1981). 
These are fully supersonic cases and are modeled using the Fabri choke method. Results are presented in 
Figure 2.4 and Figure 2.5. 


entrainment 

W2/W1 


Al/A3= 33 


♦ 



/ A1/A3-5 


/ 

♦ 


P02/P01 

Constant Area Ejector Flow. Data from Addy et al. ( 1 98 1 ). 
Conditions: Ml =2.5, T01/T02=l ., NPR=5. 


FIGURE 2.4 Comparison of Fabri choke model entrainment predictions to experimental data. 
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entrainment 

W2/W1 



P02/P01 

Variable Area Ejector Flow. Data fromAddy et al. (1981). 
Conditions: Ml =2.5, Al/A3=0.814, A2/ A3 =0.764, NPR=5. 


FIGURE 2.5 Comparison of Fabri choke model entrainment predictions to experimental data. 


The comparison with experimental data is shown to be good for a range of operating conditions. In spite 
of the fact, that these ejectors are simple slot ejectors that would not be typically used for a practical 
aeropropulsion design, they do provide confidence, that the simple control volume models are providing useful 
information. 


2.2.5 Critical back pressure computation and subsonic flow limiting conditions 

Sections 2.2.3 and 2.24 have developed two basic possible ejector flow problems, namely, back pressure 
constrained and back pressure unconstrained flows. These two flow regimes have significantly different 
modeling methods associated with them. Although the distinction that two flow problems exist is clear, at no 
point did we discuss how one chooses, “apriori”, which is the appropriate model. In this section an 
approximate model to choose the appropriate flow regime is discussed. In addition to the supersonic/subsonic 
regime transition there are a series of constraints imposed upon the subsonic operation of the ejector. 
Attempting to operate outside of these constraints will cause DREA to fail due to physical reasons. As such, it 
is of considerable interest to provide an analysis to delineate these operational regions. Though it is beyond the 
scope of this study to develop operational maps, one of the constraints is modeled and results are compared to 
DREA to show that convergence failure of DREA is indeed associated with one of the physical constraints. 

The choice of the proper condition is dependent, as one would expect, upon the back pressure or exit 
pressure level that the ejector senses. The task here is to develop a value for the exit pressure that separates the 
back pressure dependent regime from the back pressure independent regime. This value is termed as the 
critical back pressure or p cnI . The computation of this value follows. Referring to Figure 2.6, we consider the 
combined model. 
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FIGURE 2.6 Critical back pressure analysis method definitions. 


The first portion of the model involves the Fabri limited flow. The solution for this problem (which have 
been described in section 2.2.4) are the conditions associated with supersonic flow and choking in the 
secondary stream. Hence, if it were possible to estimate the exit conditions associated with this flow, then we 
would know the critical back pressure level. If the exit back pressure level is below p^, then the Fabri-choke 
is the correct solution. If, on the other hand, the back pressure level is greater than p crit then the flow is not 
properly modeled using Fabri theory and should be modeled using the back pressure dependent model. This 
criterion is summarized in Table 2.2. 


( Pe ) < ( Pcr ") 
Pot Po\ 

Back pressure independent flow 
Section 2.2.4 

( Pe )> ( Pcr ") 
Pox Po\ 

Back pressure dependent flow 
Section 2.2.3 


TABLE 2.2 Critical back pressure criteria. 


From Table 2.2, it is apparent that the key parameter is the value of p crit /poi. We may compute this value 
by reformulating the classical 1-d conservation relationships shown in section 2.2.3. These are repeated for 
convenience: 


m/5 + rrt 2 s = nu 


( 2 . 12 ) 


thlsUls + P h Als + ril2sll2s + P 2 s A2s = m3U3 + P 3 A3 


(2.13) 


rhisTois + n\2sT 02s - m 3 To 3 


(2.14) 
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and the exit pressure definition: 


P 3 = Pe (2.15) 

Combining these relationships to obtain a relationship in terms of M ls , M 2s , and the critical pressure ratio 
Pcn/Poi- Details for this reduction are presented in Appendix B. The resulting non-linear, algebraic equation 
may be written: 


O + ^T M h )7Tr (1 + y MlJ + (^-)(—)(l + Y -^~ Mb&O + 7 M 2 2s ) - 
2 Poi A, 2 


[—.<—?+ 

r-i Pot 


2jj 

y-1 


* 


[Mi 0 + ~~ Mh/i-r + Mis M 2s (—)(—)(— )2(1 + —) * 
2 P02 A, T 02 T 02 


(2.16) 


y-1 y-1 y+i 

(1 + 2—— Mi)w7) (1 + 2—— Ml)w-r) + 


Mi r— /r— /r/ + r -M- mh0pi + = o 

Pol A, 2 y-1 p 0 , A, 


Equation (2.16) relates the critical pressure ratio p C n/Poi to the aerodynamic choke quantities, Ai s , Mi s , A 2s and 
M 2s . Solution of these relationships is relatively straight forward, in that, Pcn/Poi is decoupled from the Fabri 
analyses, permitting solution of the Fabri problem and inversion (using a single variable non-linear equation 
solver) to compute the critical pressure ratio. 

At this point the limiting flow field conditions for subsonic mixing are discussed. The meaning of a 
“limiting” subsonic flow in this application describes flow fields for which the intended ejector operation; 

i.e. primary stream flow inducing secondary entrainment, thereby, satisfying the downstream pressure 
constraint fails. It is possible to identify three possible limiting subsonic limiting flow situations: 

1 . Incipient reverse flow into the secondary inlet. 

2. Choked flow in the secondary inlet. 

3. Choked flow in the exit mixing stream. 

These three cases are illustrated in Figure 2.7: 
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Potential backflow condition or 



FIGURE 2.7 Limiting subsonic flow fields. 


Our task is to determine upper and lower bounds on (Poz/PoiXA^Aj), such that there is sufficient secondary 
pressure to prevent back flow, without inducing secondary inlet choking or exit stream choking. Hence an 
estimate for (Po 2 /Poi)(A 2 /A 1 ) crit which is used to delimit various flow regimes is developed. The details of this 
procedure are described in Appendix B. 

2.2.6 Non-ideal entrainment/pumping 

The requirement that we must model non-ideal pumping effects is central to our understanding of ejector 
flows. The need for this type of modeling comes from the simplified treatment of entrainment or pumping 
employed in this project. In the previous models we introduced quasi-one-dimensional assumptions or 
isentropic restrictions. Of course real, turbulent mixing flows are strongly influenced by multi-dimensional 
effects, as well as, irreversible losses (heat transfer and viscous dissipation). In this section, reduced 
complexity closure or feedback mechanisms that model these non-ideal effects are developed. Non-ideal 
pumping effects, just as inviscid pumping modeling itself, may be conveniently divided into two problems: 
back pressure dependent pumping and back pressure independent flow. 

The closure for the back pressure dependent non-ideal pumping problem is dependent upon correctly 
modeling the multi-dimensional form of the exit pressure matching constraint. This is performed by 
generalizing the pressure constraint. It is extended by computing the integral average after the pressure has 
been corrected for viscous and multi-dimensional effects. Consider the relationship with "n" as the known 
iterative level (from the viscous mixing code output) and M n+1” as the unknown iterative level: 

P n J, =PL = \' 0 P n dA (2.17) 


Clearly, for a well mixed flow, the only non-ideal effect which will be recovered is the pressure loss effect due 
to friction and heat transfer. Fortunately, it is possible to develop an estimate of the average exit pressure 
without the burden of solving the full mixing equations for each pass. This deviation from ideal pumping is 
presented in Figure 2.8: 
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FIGURE 2.8 Back pressure dependent non-ideal flow definitions. 


Back pressure independent pumping which was modeled by ideal (isentropic, quasi-one-dimensional) 
Fabri choke theory must also be corrected for multi-dimensional effects as well as frictional and heat transfer 
effects. As in the back pressure dependent case, we estimate the amount of irreversibility and deviation from 
one-dimensional theory by using result from the viscous mixing code. These losses are presented in Figure 2.9. 


M=1 

Aerodynamic, "Fabri" choke 



FIGURE 2.9 Back pressure independent loss definitions. 


Saying this in functional form, we write: 
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(2.18) 


Now, to compute averaged pressure ratios and temperature ratios, we average using the expected relationship: 



-Vjo 1 F!’(x,y)dA 

A i S 


(2.19) 


Hence, the back pressure independent flow problem is closed, as well. Again, our goal is to develop models 
that do not require the solution of the full mixing equations and subsequent averaging. 
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The description of the algorithm above provides us with a methodology for computing the losses expected 
for the Fabri or aerodynamic choking problem, as shown in Figure 2.7. Some objections to this overall 
approach of averaging the viscous flow field must be considered: 

1. a fully 2-d viscous solution is used to correct the quasi-l-d Fabri solution. This model fidelity 
inconsistency needs to be addressed. In other words, do we need to compute an expensive, 2-d flow field 
to estimate a 1-d quantity? 

2. from preliminary experimental evidence, the effects of total pressure loss and heat transfer are expected to 
be small. Thus, the use of an expensive globally iterative algorithm may not be justified. 

With these objections in mind, the need for a simplified method for the estimation of the loss becomes 
apparent. 


2.2.6. 1 Back pressure dependent regime loss/heat transfer model 

As described in section 2.2.3, pumping for the subsonic problem is modeled by demanding conformance 
to the 1 -d/average exit static pressure condition: 


= \‘ 0 p n dA 


F = Ple-P, 


exit,amb. 


( 2 . 20 ) 


Since the model already iteratively matches p„it=Pambi e nt assuming complete mixing, what remains to be 
computed, is a correction for the profile. This is especially important for mixing processes where significant 
regions of supersonic flow may exist: 

In this type of problem or in a shock down problem, the lack of mixing will strongly influence the average 
pressure field and thereby the entrainment. In contrast, the entrainment for a fully subsonic ejector will be 
significantly less sensitive to the degree of mixing. As such, we will always model fully subsonic ejectors using 
the ideal mixing approximation. 

Returning to flows with significant regions of sub sonic flow, i.e. M t >l, we consider how the exit static 
pressure profile may be estimated. This could be done by integrating the profiles provided by the full mixing 
equations, but as pointed out, this is computationally prohibitively expensive. Consider alternatively the 
development of an approximate differential equation written for the average pressure. We start by defining the 
dimensional relationship (essentially a finite difference approximation): 


dPave Pave , 3 Pa 

dx L 


Pave = JoPd4 


char, mix 


( 2 . 21 ) 


where L char miK is a length characteristic of the mixing and p lv0 is the ideally fully mixed value. Clearly for a 
very long ejector, the average pressure must asymptotically approach this value. As is shown in section 2.4., it 
is reasonable to associate this length scale with: 


-'char. mix 


V g . ‘V I 


U m H ! 


( 2 . 22 ) 


Non-dimensionalizing and substituting into the previous equation yields our final differential equation and the 
initial condition: 
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= ax(p avei -p are ) 


Povc(°) = PA + 0 - 


(2.23) 


JPave 


dx 


where a*x is the dimensionless effective turbulent viscosity. This is a linear differential equation which is easily 
solved to yield: 


= j + [PaJVl _ j] exp(-^l) Pa J0) = pA+( 1 " K )Pz (2-24) 

Pave t 3 Pave ,3 


With this solution, it is possible to estimate the effect of mixedness upon the secondary properties. Since now 
the exit constraint takes the form: 


P~P ove (i shroud') ~ P exit ,amb ~P ave ^ L shroud )J Pexit.amb 

P ave. 3 


(2.25) 


Equation (2.25) drives towards zero at convergence. 

Entrainment or pumping is a strong function of the shroud length. A poorly mixed problem will not satisfy 
the 1-d exit constraint, and will fail to give a reasonable estimate of pumping. We note that pumping will still 
occur, by a viscous, jet like entrainment, but this limiting effect is not modeled here. This is not a critical 
limitation, since any practically designed ejector will seek to maximize mixing, typically by using a complex, 
lobed mixer. 


2.2.6.2 Back pressure independent regime loss/heat transfer model 

A model based upon extensions of the inviscid analyses developed previously for back pressure 
independent flows (to include loss terms) and approximations of the more complete mixing model are 
presented here. The solution procedure may be stated as: 

• estimate the choke location 

• compute the asymptotic or fully mixed losses, which will be needed for the loss differential equation 

• estimate the turbulent viscosity, which is also fundamental to the loss equation 

• solve the loss differential equations 

• substitute the loss terms into the extended Fabri choke equations. 

Since this is a rather complex problem, a useful way to visualize this progression is to write the functional 
forms of these equations with their variable dependencies. Further, due to the simplifying assumptions applied 
to the associated differential equations, this model is completely algebraic. As such, we may solve it with any 
one of our non-linear inversion algorithms. See Section 3 for a discussion of nonlinear equation solution 
algorithms. Consider then the vector equation written in functional terms: 
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X c (M v M u Jou>T^H s ) 

To.. = T 0< MvP m ,T 0 „T, 2 ,M x ,M 2 A x ,A 2 ) 

= hf 3 (P 0 i,P 02 ,T 0X ,T 02 ,Af t ,M 2 'A x ,A 2 ') 

Po.*>=PoA P Ol’ P 02’ToiJo2>M l ,M 2 A’ A 2>M 3 ) 

a =a'(P 0 l ,P 02 ,T 0 ] ,T 02 ,M x ,M 2 A ] ,A 2 ) 

^Ols ~ Pois(Po\’Po2’Po,«>’ a >* C ) 

Po2s = Po2s(Poi’Po 2 ’Po,n’ a »-^ c ) 
Tou^n u (T nt T nt T 0 ^a\X e ) 
To 2 s=T Q2s (T (n ,T 02 ,T 0 ' X ,a',X c ) 

■ (P 01 , P 02 ,P 0U , P 02s , To, , T 02 , T ou , T 02 s ,M x , M 2 ,M u ,M 2s ,A x , A 2 ) 

\momentum (P 01 , ^02 5^01j>^02j5 ^0P^02’ ^01* ? ^02 j 5 M \ , M 2 , A/ l5 , M 2s , 4 , A 2 )J 


mass 


= 0 (2.26) 


and the (1 1) unknowns: 


(*C 




^ fl 


1 01 s 


‘02 5 


7" 

I ou 


l 02s 


K m 2 )’ 


(2.27) 


Solution of these equations is simplified by the fact that the first nine equations in the system may be written 
explicitly. Thus, although this system is formally an eleven equation system, it may be solved using a two 
variable inversion. Broyden’s method is used here (Burden and Faires (1993)). Also see Section 3. 

The “workhorse” analysis in the above system, are the mixing equations for the total pressures and 
temperatures. All the analyses have been taken from previous sections (either with some simplification or with 
extensions). Since these governing ordinary differential equations are new, their development is outlined here. 
This starts with the “model” equations for mixing of any scalar quantity. The notation, <f> is chosen to 
emphasize the analogy with the main governing equations: 


dfa_ _foc , 

^ ave Tj(rl r«>) 

ax H 


(2.28) 


u b* (A 


(2.29) 


where “kx” represents v eff the turbulence model and (j)* is the asymptotic fully mixed value. The subscripts, of 
course, refer to the appropriate stream. Non-dimensionalizing yields: 


d ± i 

dx 


= ax{<p i -</>J 


(2.30) 
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= ax(<j> 2 -<j>J 
ax 


(2.31) 


Since these equations are linear, they may be immediately solved. It is interesting, though, to subtract these 
equations to yield: 


zhl 

dx 


= -a'x(<f> i -<f> 2 ) 


(2.32) 


In this form the mixing form of these equations is more apparent. 

Solving equations (2.32) and (2.33) provides estimates of the total conditions within a stream at any place 
in the flow field. Consider the general solution: 




(2.33) 


Thus we have access to the total quantities (P 0ls , Po 2 s* T 0 i s , T 02s ) at any point in the flow field. Although this 
methodology is rather simple, it does contain both an approximation to the turbulent mixing, as well as, the loss 
mechanism through the asymptotic value. 

The Fabri governing equations must be extended to include the loss terms as computed above. This is a 
straightforward derivation where loss terms are now retained. Consider the extended equation for mass 
conservation: 
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and the global momentum equation: 


(2.34) 
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(2.35) 


The extension terms are apparent in these equations as well as the closures (P 0 i s , P 0 2 s> T 0 i s , T 02s ) which were 
described previously. 

The above paragraphs provide a description of the algorithm that is to be used to model the form of non- 
ideal mixing effects on the pumping. This provides us information on the flow rate and hence, other static 
quantities. Unfortunately, the above equations do not give us any idea of the convergence rate of the method or 
conditions under which convergence (or divergence) might be expected. Of course, we can make some 
estimates of these effects by referring to general fixed point iterative methods ((Johnson and Riess (1982) or 
Gerald (1980)). Some areas that might be considered include: 

1 . General fixed point convergence rate linear 

2. General fixed point convergence rate criterion 

3. methods to improve convergence rates (such as Aitkens method (Johnson and Riess (1982)) etc. 

Fortunately, since the model was described by the reduced set of equations, the convergence rate is the same as 
Broyden’s method, which is superlinear. 


2.3 Inviscid flow field analysis 

As indicated, the viscous mixing region solver is a steady state, single pass marching algorithm. Although 
the basic mixing process is formulated in terms of conservative flux quantities (mass flux, momentum flux and 
energy flux) which are Mach number invariant, the local primitive variables, such as, velocity, temperature, 
static and total pressure, density and Mach number are certainly subsonic/supersonic dependent. Since the 
supersonic or back pressure independent ejector problem, is by its veiy nature, a transonic flow; delineation of 
regions of supersonic or subsonic flow must be provided for the mixing code primitive variable inversion 
analysis. Estimates of the flow regime are supplied using a set of extended, classical inviscid, 
supersonic/subsonic slipline interaction algorithms. These are described in this subsection. 


2.3.1 Slipline interface analysis 

Considering the most general problem, we note that a supersonic gas jet issuing into a co-flowing medium 
has an almost periodic structure for small stream pressure variations. In this section, formulas that describe the 
periodicity and wavelength of the slipline are developed for two-dimensional and axi-symmetric, supersonic- 
subsonic jets. Classical inlet and streamline perturbation solution methodologies are limited to small stream 
pressure imbalances. To increase the range of validity of the streamline perturbation solution, a modified 
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strained coordinate technique is introduced that shows good comparison with experimental data over the Mach 
number and pressure ratio p,/p 2 range of 1. 2-4.0 and 1.0-20.0, respectively. This more general relationship 
provides a limiting form of a supersonic ejector flow field, since the location of the first maximum expansion 
corresponds to the location of the aero-dynamic or Fabri-choke analyzed in the previous section. 

An early class of models for free jet flows were developed by Prandtl (1904) and later by Pai (1952) and 
(1954) based upon perturbations about the inlet quantities, i.e. u=U,+Ui’. Although straightforward to use, these 
models have restrictive assumptions and yield relatively poor results. An improvement was made by Pack 
(1950) who used a perturbation about the velocity field at the slipline. This model shows significantly better 
comparison with experimental data but is still only valid for small pressure ratios between the jet and ambient 
flows. 

In this section, the fundamental equations are merely quoted. The derivations are deferred to Appendix C. 
The appendix starts by briefly describing the models of Pai and Pack for reference use. Then, a modification 
to the streamline perturbation model is derived using a strained coordinate technique that extends the range of 
application to larger pressure ratios for jet flows. A comparison is performed to evaluate the new model for 
available large pressure ratio free jet experiments. Having helped establish the validity of the relationships for 
free jets, their use for the ejector flows is developed. 

The principle parameter is discussed is x c , which is the location of the first minimum in the 
expansion/compression system resulting in the flow. This value for both axi-symmetric and two dimensional 
flows is described by the equations below: 

•axi-symmetric: 


~Z = 1.22P eff -.22P ts - 1.22(1 + -$ m )P h -.22P ls 


(2.36) 


•two-dimensional: 


~zr = 2 P e jj - P ls = 2(l + ~: ZJ)P ls - P ls 


(2.37) 


where p ls =(M, s 2 -l) 1 ' 2 and ^=P, S 2 [1-U,/U, S ]. Equations (2.38) and (2.39) are the final modified solutions for the 
critical slipline location. These relationships with the lowest order terms based upon control volume theory 
comprise our solution for the location of the first minimum of the slipline. By symmetry, half this distance is 
the local maximum of the primary expansion and corresponds to the location of the Fabri-choke. Results are 
presented in Figure 2.10. Further comparisons are presented in De Chant (1996). 
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FIGURE 2.10 Comparison of new strained coordinate analysis interface geometry (first critical) location, Xc 
with inlet perturbation and experimental correlation of Love (1959). 


A simple linearized, potential flow model for unconfined/free jet, axi-symmetric and two-dimensional 
flows has been developed based upon a strained coordinate extension to Pack's (1905) streamline or slipline 
perturbation. 

This analysis provides analytically based estimates of the slipline geometry for co-flowing supersonic and 
subsonic streams. The slipline is described by the wavelength or the streamwise location of the first minimum 
of slipline displacement. Free jet data was used to test this model and compare with other classical 
relationships. 


2.3.2 Confined flow modifications 

» 

The prediction of the critical location or wavelength for unconfined flows developed here provides a 
limiting case estimate for ejector flows which are important for aerodynamic and industrial applications. 
Additionally, this model provides information about the efflux of a supersonic jet, which is important for 
aerodynamic noise generation problems. Of course, though, ejector/mixer nozzles are by their very nature 
confined or internal flow problems. As such, modifications to the basic inviscid flow analyses are necessary. 
This modification takes the form of estimating the base flow quantities in a different manner. Consider the 
expansion 


u dd> 

— = ; + 

U is dx 


// U 1 \ 


(2.38) 
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where, U )s , is the velocity at the interface. Note the local definition of <|>, which here represents a small 
disturbance velocity potential. 

As explained in Appendix C, the slipline velocity U )s , was estimated using a simple static pressure 
matching methodology for the free jet problem. To account for the internal flow effects, though, it is necessary 
to apply a control volume methodology. Consider the two fundamental problems for an internal flow problem: 

(1) Large area ratio, unknowns: M ls , M 2s , U ls 

(2) Small area ratio (Fabri problem), unknowns: M ls , M 2 , U| S 

The necessity for two problem closures needs to be explained. The conditions for a free jet, are based solely 
upon local static pressure matching. Consequently, problems for A 2 /A, »1 are analyzed using static pressure 
and mass conservation. Unfortunately, the local 1-d control volume assumptions are poor approximations for 
sufficiently small area ratio’s, A 2 /A h . Therefore, this approximation cannot be expected to apply for A 2 /Ai ~ 
1, and hence provide an alternative formulation. 

Starting with mass conservation which is common to both problems, we write: 
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(2.39) 


and static pressure matching which is used for large area ratio problems A 2 /A t »l 
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(2.40) 


or the global momentum equation used for, A 2 /A 1 «l. 
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These relationships provide closure for the internal flow Mach numbers and area ratios. 

Since the flow is locally assumed to be isentropic, it is possible to compute the slipline velocity: 


U i Mi 

U is Mis 

Contained as a degenerate case is the free jet solution: 
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l + ~-M] 


(2.42) 


Mi = iTl^ 1 + L r r - 1] (2-43) 

y-l 2 p 2 

With these low order forms available the analysis proceeds in a straightforward manner, namely solving and 
the Fabri choke location using equations (2.36) or (2.37). 
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2.4 Parabolic mixing /low analysis 


The second major portion of the flow chart, Figure 2.1 considered is the turbulent compressible mixing 
analysis. In this section the character and solution of the mixing governing relationships which were developed 
in Section 1 are discussed. 


2.4.1 Conservation form mixing differential equations 

A reduced set of partial differential equations was derived in Section 1, i.e. equations (1.19), (1.20) and 
(1.22). They are a decoupled set of parabolic, scalar, heat equations. Equations (1.19), (1.20) and (1.22) may 
be placed in a simpler, “canonical” form. Next, numerical and classical analytical solution methods are 
developed and the inherent limitations of both methods are discussed. Finally an alternative to these solely 
numerical or analytical solution methods which combines the best features of both numerical and analytical 
methods is described. 

2.4.2 Solution of mixing partial differential equations 

From the previous section, the governing equations may be approximately written for the conservation 
quantities in terms of the general linear parabolic equation: 


dtp , .d 2 $ ♦ d 2 <f> 

— = a(x) — rr = a x—T 
dx dy 2 dy 2 


(2.44) 


with the boundary conditions: 


d<f>(x,0) d(j>(x,l) 

dy dy 


(2.45) 


with the initial condition: 


<(>(0,y) = 


$1 to 

<t>20 


0 < y < 
h s <y<f 


(2.46) 


and <|> is defined by: 


(pu 2 + p) 


<fi(x>y) = 


puH 


\ pu ) 


(2.47) 


The complex turbulence model (v e n/U) 0 G/e has been replaced by some function a(x). This functional closure 
will be described in detail in Section 2.5. Additionally, the superscript and the “0” subscript denoting a 
dimensionless quantity and level within the perturbation expansion, respectively, have been suppressed for 
simplicity. 
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For our turbulence model it will appropriate to make an additional local assumption, that the function a(x) 
may be approximated by the linear function a(x)=a’ x. This will give us an analogous formulation to that of 
Reichardt. See Schlichting (1979). The choice of the locally linear functional form is justified by our 
development of the near field turbulence model. Again, this assumption will be justified in Section 2.5. 


wall shroud 


secondary 


primary 


1 


centerline 


FIGURE 2. 1 1 Mixing analysis singularity definitions. 


The above governing equation with associated boundary and initial conditions are on the surface a 
relatively trivial set of equations. As we shall see, though, the discontinuous nature of the initial condition, as 
shown in Figure 2.1 1, will cause considerable difficulties both for numerical and analytical solutions. With this 
in mind, it will be of interest to “naively” solve the governing equation. It is linear so one may immediately 
write an analytical solution. The straightforward eigenfunction expansion (Haberman (1983)) solution is: 

<P(x, y) = <P + ^(^kj- 4>2q) 2-T sin f nnh s ) ZQ^n7iy)e a(nK), l (2.48) 


where: 


<t>= < t > i 0 hs + (1 -h s )<j> 20 = 1 


(2.49) 


This equation is formally an exact solution to equations (2.46)-(2.50). Unfortunately, this equation is very slow 
to converge for x«l, indeed truncation error is on the order of 0(1 /n), and the solution suffers from Gibb’s 
jump departures (Haberman (1983) and Weinberger (1965)). Unfortunately, x«l is precisely where good 
converge is necessary since the concentration in this project is on short shroud ejectors. 

Since the eigenfunction expansion (Fourier series) method is both inefficient and yields very poor results, 
alternative methods must be considered. The most obvious solution from the standpoint of efficiency is to 
apply a strictly numerical solution. Early in this work, a high accuracy (0(Ay 4 ), 0(Ax 2 )) numerical solver 
based upon a variable step, compact differencing scheme (Kreiss differencing) in the cross-stream and variable 
step, fully implicit differencing in the streamwise direction, was applied to this problem. This model is 
described in detail in Section 3. Although, the solver which achieves 4th order accuracy on three-point 
support by simultaneously computing the function, and two derivatives, was of high accuracy, it could not in 
general, resolve the discontinuous nature of the input profile. This limitation was shown clearly when the 
divergence theorem constraint, i.e. mass conservation which demands: 
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\‘ 0 <p(x,y)dy = 1 


(2.50) 


Unfortunately, this constraint was not properly honored by the strictly numerical solver due to the singular 
behavior near the step input. Clearly a methodology which is not even conservative is unacceptable. 

The poor performance of the strictly numerical model may be analyzed further. Though Section 3 
describes the compact (spline-on-spline method) in detail, it is useful to discuss the truncation error briefly 
here. Considering the implicit relationship governing the second derivatives we find the error terms: 


error =-j-(A J+ iAj (A J+ i - Aj))<f> W + 0(A 4 )<!> {b) 


(2.51) 


which is formally fourth order for constant grid spacing. Unfortunately for very rapid changes in <j>, the high 
order derivatives, such as, <}> (6) are very large. In fact they may be as large as 0(1/A 4 ) which would yield a very 
poor solution with errors 0(1)! Indeed, this was observed. 

Clearly, naive analytical solution or numerical solution of this problem will lead to poor performance. 
What is requires is an analysis capable of dealing with extremely rapidly changing functions. Fortunately a 
solution method from classical analysis is available which is based upon distribution theory rather than 
continuous functions. This is the method of Green’s function expansions. We may actually, more usefully 
derive the solution using a cosine transform, i.e. even function, semi-infinite Fourier transform and the method 
of images (Haberman (1983), Churchill and Brown (1978)). The result, which is merely quoted here, is derived 
in Appendix D. 
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(2.52) 


Note, that for x«l this relationship recovers the step input. The solution could have also been derived directly 
using infinite space Green’s functions and the method of images (Haberman (1983)). 

Although the above relationship is exact, and does not suffer from the near field (small "x") limitations that 
the eigenfunction expansion solution did, it is still in the form of an infinite series. It may be noted that, in the 
near field, the solution converges extremely quickly. This is desirable since most of our interest is in the near 
field. Although short shroud ejectors are of interest in this study, it is necessary to be able model far field 
behavior using a reasonable number of terms in the series. A method, which invokes the best features of the 
analytical near field expansion and an efficient finite difference approximation, is developed subsequently. 

It is important to quantitatively estimate the number of terms required to satisfy the governing equation to 
some pre-determined level. This knowledge will make it possible to balance our numerical effort against the 
analytical effort, since a combined problem is ultimately to be solved. This may be done most conveniently by 
considering the derivative of the flux at the wall, i.e. 5<j>/5y. In particular, it will be useful to consider the size of 
a single term in the Green’s function expansion. The choice of using the boundary condition as an indicator of 
convergence, is especially appropriate, in that the Green’s function expansion is very accurate near the step 
function (Haberman (1983)), with reduced accuracy far from the singularity. Hence a boundary condition is a 
reasonable requirement used to estimate convergence rate. The single term is written: 
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(2.53) 
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Selecting y=l and concentrating on the rate of exponential decay (for the first term in the above equation), and 
solving for n yields: 
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(2.54) 


This gives us an estimate of the required to obtained a pre-specified error for dtyldy. Two trends may be 
observed: 

1 . As x or a* increases, the required number of series terms, “n” required will increase and 

2. As the difference in initial condition <})jo“<}> 2 o increases, the required number of terms will also increase. 

Numerical experimentation, shows that a maximum number (even for very long ducts) of terms in the Green’s 
function expansion is always less than six. This is acceptable, since the problem is solved using a hybrid 
numerical/analytical method, where the series truncation error is corrected by invoking a numerical/analytical 
decomposed problem. Again, hybrid method is described subsequently. 

A combined numerical and analytical solution of the problem is easily effected, since the governing 
equation (2.47) is linear. Consider the decomposition: 


<K*> y) = 'Pan ( X > y) + Pnum ( X > 3 ^ 


(2.55) 


with the governing equation: 


dftnum _ n / Y \ d $num _ * y & num 

dx W dy ’ dy 1 


(2.56) 


with the new boundary conditions: 


&num <K X >0) _ d an <p(X,0) 
dy dy 


d num (x,l) _ danftxj) 
dy dy 


(2.57) 


and the new initial condition: 


<PnuJ0,y) = 0 


(2.58) 
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The overall solution is then, of course the sum of the two functions. Notice, that the new numerical initial 
condition equation (2.58) is essentially trivially, and very easily represented numerically. 

Before concluding the mixing analysis, it is worth noting, that the singularity that is shown in Figure 2.1 1 
is caused by our mathematical representation of the near field conditions. Physical, a discontinuous initial 
profile cannot exits in the presence of physical viscosity. In reality, the flow has wall boundary layers and a 
“blunt” tip zone (and therefore, a locally subsonic elliptic region) very near the splitter plate as shown in Figure 
2 . 12 . 



FIGURE 2.12 Details of the physical flow field near the splitter plate. 


We chose not to attempt to model these effects for the following reasons: 

1 . Inlet conditions are complex and difficult to model at the level of fidelity which the mixing model has been 
developed for. Additionally more information than will ever be known at the initial design level would be 
required to properly model the inlet flow. Figure 2.12 is clearly reminiscent of classical “triple deck” 
theory (Van Dyke (1975)), with all the attendant theoretical and analytical complexities associated with 
that work. 

2. It is our intention to connect this mixing model to thermodynamic design codes, for which profile 
information is not known, whereas 1 -d inlet information certainly is. 

3. Modeling of the near field would require an expensive elliptic methodology, completely inappropriate for 
the model that is developed in this study. 

For these reasons, the decision has been made to model a virtually (and mathematically exactly) an input 
step discontinuity. This is effectively a Heaviside step input (Haberman (1983)). The purpose of this 
paragraph, though, is to state that the inlet discontinuity is a modeling feature. Colloquially stated, physical 
fluid flows are “smart enough” to prevent the existence of singular behavior, we have introduced this 
idealization for the reasons stated above. 

2.4.3 Primitive variable decoding algorithm 


In this section, the inversion of conservative flux quantities, <p, to “primitive” variable such as velocity, 
pressure, temperature, density, Mach number etc. is developed. 
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2.4.3. 1 Complete inversion 


The solver described in the previous section is written strictly in terms of the conservative flux quantities: 


<f>(x,y) 


pu +p 
puH 
P* . 


(2.63) 


To convert, these fluxes into primitive variables, such as (M, u, p, T, p), a local one-dimensional 
approximation is applied combined with the definitions of the fluxes themselves to compute the primitive 
variables. Consider, for example, the velocity may be recovered from the flux values using: 

F ( u ) = [ pu(x, y)]( u f - [( pu 2 + p)(x, y)]u + ^~[ puH(x, y) ] = 0 (2.64) 

zy Y 

A somewhat more convenient form for analytical purposes of equation (72), which is in terms of the 
conservative fluxes and the Mach number, may be written: 
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(2.65) 


These equations contain a considerable amount of physics. The two roots of the quadratic equation (in terms of 
M 2 ) correspond to supersonic and subsonic solutions for the ejector flow field. Hence, at every point in the 
flow field, the two basic ejector solutions are contained. One can also show that where the flow is choked, 
M=l, implies a single solution. This corresponds to a zero discriminant in equation (2.64) or (2.65). The 
choice of the appropriate solution is determined by inviscid flow field solver, which locates the Fabri-choke 
and the subsonic/supersonic interface. The very close correspondence to shock capturing methods used for 
Inviscid Euler type solvers (Anderson et al. (1984)) should also be noted. With this similarity in mind, it is 
apparent, that this equation has strong normal shock solutions embedded within it. This is consistent with our 
ejector analysis, which for negligible secondary flow must recover the classical-one-dimensional normal shock 
relationships. 

Although equation (2.65) provides a method for recovering primitive quantities from the conservative 
fluxes, for the steady, viscous problem with embedded subsonic regions, this formulation is not without 
problems. This problem becomes clear if one considers Figure 2.13: 
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FIGURE 2.13 Characteristics of complete and modified primitive variable inversion relationships. 


From Figure 2.13, the double valued nature of the function C 0 is apparent. Recall that C 0 is the 
dimensionless ratio used in equation (2.65). As the graph implies, it can also be written solely as a function of 
Mach number: 
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which is merely equation (2.65) rewritten. This double valued nature, although physically significant in the 
streamwise direction (normal shock down shock or Fabri expansion) causes considerable trouble in the cross- 
stream, or “y” direction for mixed supersonic/subsonic flow fields. 

Consider the mixing of low speed (large C 0 ) subsonic fluid (subsonic branch) fluid with supersonic fluid. 
Mixing causes a decrease in C 0 in the subsonic region which, correctly, decodes as an increasing Mach number. 
Unfortunately, this large value of C 0 will also increase the C 0 within the supersonic portion of the mixing zone 
(the flux quantities are monatomic and hence the value of C 0 is also). This incorrectly decodes as an increasing 
Mach number on the supersonic branch! This is represented in Figure 2.14. 

This cannot be correct since the Mach number must monatomically decrease from supersonic to subsonic. 
Other primitive quantities are similarly effected, including local static pressure spikes. 
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What has happened? The “failure” of the model is caused by the singular behavior of the inversion 
function for M«l, or incompressible flow. Equation (2.66) is not bounded for M-»Q, which in turn is a 
remnant of the fact that the state equation and the incompressible assumption 

p - pRT p = const. (2.67) 


cannot simultaneously hold. Of course, for a single regime, i.e. fully subsonic or fully supersonic flow, the 
singularity is irrelevant except for M=0, exactly. Hirsch (1988) and Anderson et al. (1984) also discuss a 
related issue, though, this is more of a question of computational expense when applying compressible flow 
algorithms to virtually incompressible problems. 


large (subsonic) Co mixes with smaller 
(supersonic) CO. This decodes correctly 
in the subsonic region (branch), but decodes 
incorrectly in supersonic region (branch). 


large value of Co 



large Co, 
susonic brach 


smaller Co, 
supersonic branch 


1 

1 

Mach number 


FIGURE 2.14 Primitive variable inversion failure for mixing transonic flow fields. 


The problem encountered here is closely related to the departure solution phenomenon encountered for 
mixed subsonic and supersonic flow fields when solved by single pass parabolic Navier-Stokes solvers 
(Anderson et al. (1984)). In this type of primitive variable formulation, elliptic influences of the subsonic 
pressure field cause the development of negative eigenvalues and hence poorly posed, and unphysical 
solutions. These take the form of unphysical pressure fields. It was hoped, by the author, that this problem 
would have been eliminated in this study, by the choice of conservative variables, since the scalar eigenvalues 
of (2.44) a simple heat equation are indeed positive, but at a cost of a non-physical inversion. 
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2.4. 3.2. A model problem for primitive variable inversion 

Since the inversion problem is so closely related to the ill-posed behavior exhibited by space marching 
version of the parabolized Navier-Stokes equations (Anderson et al. (1984) or Vigneron et al. (1978)), it is 
useful to explore the use of transformations and subsequent inversion further. An appropriate way to do this is 
to consider a model equation (which is closely related to our problem). Consider the equation: 



dM 

dx, 


v_cff c?M_ 
U dy 1 


( 2 . 68 ) 


This is a scalar analogue to equations (2.44)-(2.47). Equation (2.68) exhibits two regimes a supersonic 
regime (M>1) and a subsonic regime (M<1). Examination of the lead coefficients indicate that equation (2.68) 
will not yield a physical relationship for the subsonic region, since this would imply, that initial conditions 
would grow in an unbounded fashion as the solution marches forward in space, i.e. a negative eigenvalue. 
Alternatively stated, for subsonic flow with local linearization the “viscosity” 
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(2.69) 


would be negative! This is contrary to our physical understanding of diffusive phenomena. Looking at a flux 
vector in this case “Mach flux” in the cross stream direction: 




V dM_ 
U dy 


(2.70) 


would lead to the conclusion, that Mach number increases from low potential to high. . This would violate the 
second law of thermodynamics. 

Since it is apparent, that the current formulation leads to unacceptable behavior for subsonic flow, it is 
necessary to consider alternatives. The alternative developed in this study was the recasting of equation (2.68) 
into conservation form, and the approximate modification of the viscous term. This lead to: 


3 _ 

dx 
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U dy 2 



(2.71) 


In this form, the governing equation is unquestionably well posed at least with respect to the linear diffusion of 
M+l/M=G mode i. This modification is closely analogous to our linearization as shown by equation (2.44). 
Indeed by using the diffusion (or mixing) equation to predict G^i at all points in the flow field, a model 
inversion polynomial (analogous to equation (2.47)) may be formed: 
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G mode/ = M + 
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^mod el ^ [<J 2 mode/ — 4^ 


(2.72) 


which has the desired two roots corresponding to subsonic and supersonic shock phenomena. Unfortunately, 
though, well posed in terms of the linear diffusion operator and retaining shock phenomena does not imply well 
posed in a physical sense. This failure of the transformation becomes apparent when one considers a plot of 
Gmodei with respect to M. Once again this is analogous to Figure 2.13. 



Mach Number 


FIGURE 2.15 Characteristic inversion relationship for model problem. 


From Figure 2.15 and motivated by our previous discussion, the limitation inherent with change of 
variables may be discussed. A special, but insightful problem is illustrated in the Figure 2.15. Here a subsonic 
Mach field left branch is mixed with a supersonic Mach field, right branch. If the case, where G modc \ is exactly 
the same, denote by the dashed line, is chosen we have the situation where Gm^i is constant, i.e. no mixing and 
yet in the untransformed problem we have a subsonic/supersonic mixing problem! Clearly, the transformation 
has failed to honor the monatomic behavior of the mixing problem. This is because the nonlinear 
transformation is double valued. It is apparent, that a double valued function will not satisfy the physical 
requirements needed for a monatomic mixing problem. This can be reinforced if one considers the G m M flux 
vector: 
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So, in summary, the transformation yields a well posed diffusive problem in terms of G mode i, but this does 
not guarantee a well posed formulation in terms of the primitive variables. Clearly, a non-linear transformation 
does not eliminate the difficult problem of modeling viscous, transonic flows using space marching methods. 
In passing it may be noted, that a potential solution to this dilemma is the introduction of an unsteady term, 
effectively modifying the form of the governing equation: 
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(2.74) 


Equation (2.74) is effectively, a hyperbolic-parabolic equation. It would require multiple pass relaxation to 
drive to a steady state; a highly inefficient procedure relative to single pass space marching. It is not clear, that 
even this formulation will lead to a well posed problem. 

The basis for this concern comes from a simple analysis of this problem. One may linearize, equation 
(2.74), and solve by using separation of variables and the method of characteristics. The most important result 
of this analysis, is that for a supersonic problem the characteristic of the hyperbolic portion will take the form: 

* = [ 1 - 7 ^ 1 ' +c ° (275> 

M M 


This is a right running wave in the x-t plane. As per our intuition, signals (or initial condition information) are 
sent from upstream to down stream. Thus, this problem is considered to be well posed. 

Alternatively, considering the subsonic problem, we must write: 

X = -[1--!U + C 0 [l--^-]S0 (2.76) 

M M 


which is viewed as a left running wave in the x-t plane. This type of propagation does not follow our intuition, 
in that it demands that the initial condition information be propagated downstream to upstream. Hence, under 
these conditions the inlet conditions (initial conditions) propagate right out of the solution domain. This is a 
consequence of the elliptic-like behavior (for a true elliptic problem, one would expect to have imaginary 
characteristics) where downstream influences propagate back upstream. This limitation is a direct consequence 
of treating the pressure term as a part of the flow field for subsonic regions. See Lighthill (1953) for a related 
discussion. 

The solution to this problem used for parabolized Navier-Stokes (PNS) flows is to approximately, restrict 
the size of the static pressure gradient within the subsonic region. (For the model equation, this requires 
neglecting the 1/M term within the subsonic region of the flow. Fortunately this is a much better physical 
approximation for the real problem) This has the effect of modifying the non-linear function to be a single 
valued relationship. One the most successful methods was pioneered by (Vigneron et al. (1978)). By analogy 
with this work, a variation of this method has been applied to this analysis to define a new inversion technique 
that is viable for mixed flows. This method is developed subsequently. 
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2.4. 3. 3 Modified inversion 


The modified inversion method requires limiting the size of the static pressure term in the inversion 
method of equations (2.65) and (2.66), within the subsonic region only. Hence, the momentum flux is 
modified: 


G mod =P" 2 +COp 


yM 2 

( 0 = r 

\ + (y-\)M 2 


(2.77) 


where, ©, is the “Vigneron” parameter (Vigneron et al. (1978)). The form of the parameter, co, applicable for 
simple mixing flows is derived in Appendix E. Forming the analogous non-dimensional parameter to equation 
(2.66) yields: 
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(2.78) 


and the corresponding inversion polynomial: 
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(2.79) 


Plotting equation (2.79) in the Figure 2.13, repeated here as Figure 2.16: 
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FIGURE 2.16 Characteristics of complete and modified primitive variable inversion relationships. 


and observation of the form of the equation indicates that the this dimensionless parameter and the resulting 
inversion relationship are bounded for M->0. Further, the inversion is single valued and will prevent the 
unphysical behavior exhibit previously. This type of formulation will be used to model mixing flows. For this 
parameter to be available throughout the flow field requires that a new differential equation in terms of the 
modified momentum flux also be satisfied. Hence for mixed flows the vector of conservation quantities is 
now: 




' pu 1 +p' 
puH 
pu 

\pu 2 + cop) 


(2.80) 


For ejector flows, where the subsonic stream must accelerate, a slight modification of this model is 
necessary. This is due to the fact that the subsonic stream must be able to accelerate (due to the Fabri or 
aerodynamic choking phenomenon) to supersonic speed. Within the scope of one-d flow and the single valued 
inversion that has been developed, it is not possible to achieve this acceleration. The approximation introduced 
by the modification cannot be valid beyond the sonic point. The momentum field approximation, equation 
(2.77), must be relaxed and the complete momentum field is applied, for supersonic regions of the flow field. 
This is possible, since the governing equation for the complete momentum, pu 2 +p, is also solved. Hence it is 
“blended” into the flow field. This blending is accomplished using the function: 
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where G=pu 2 +p and G c n=pu 2 +cop. This function is constructed to smoothly introduce the full momentum G. 
Indeed for x»x c /2 i the choke location; G is completely recovered. Although the function is empirical, the 
basis for it is consistent with the first order approximations needed to model transonic flow using steady, 
marching equations. 

Another closure that may be used at the Fabri choke point is to merely “jump” the solution to supersonic 
and complete inversion. This has been done successfully for several test cases and is described in Section 5. 
The inherent limitation in the governing equations, i.e. that wave phenomenon such as expansion and 
compression wave are not modeled is discussed in Section 5 also. 

2.4.3.4 Flow regime and thermodynamic restrictions 

Due to the existence of two possible solutions corresponding to supersonic and subsonic flows, 
determination of the locally appropriate solution type is critical. Actually the use of the modified pressure 
field, somewhat circumvents this problem, in that the modified inversion is single valued. Even with this 
modification, though, it is necessary to determine which flow field regions will be considered supersonic, since 
we always use the unmodified inversion method for supersonic regions (just like the classical space marching 
PNS problem) and use the modified pressure field only in the subsonic region. As such, the first situation: a 
supersonic and subsonic stream interacting due to turbulent mixing is considered. 

Determination of flow field type for the mixed problem proceeds by considering the required value for the 
dimensionless parameters, C 0 and C 0imo d* Due to the formulation chosen, at the transonic point, these two 
relationships have the same value. This in fact, gives us the critical value for C 0 and C 0t mod* 


C 0 (M=l) = C 0mod (M = l) = 


2(r 2 -i) 

2 


(2.82) 


This intersection point is shown clearly in Figure 2.16. Equation (2.82) is utilized by choosing a subsonic 
solution for (and the modified inversion polynomial) C 0 , mo d ^ C 0 , mo d(M=l). Conversely, one must choose the 
supersonic solution (and full inversion polynomial), for C 0>mo d ^ C 0 ,mod(M=l). 

In addition to mixing of two streams, it is necessary to be able to identify thermodynamic restrictions 
placed upon both fully supersonic flow and fully subsonic flows. For these types of flow fields, the complete 
inversion polynomial is used. To understand the limitations on Co we may refer to Figure 2.16 or examine the 
solution types possible for equation (2.65), repeated here: 
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This second degree polynomial, in terms of M 2 , has the classical restrictions on the number of roots that are 
solutions which can be ascertained by applying Descartes Rule of Signs (Johnson and Riess (1982)). Here, we 
are looking for sign changes in the polynomial coefficients, indicating the existence of the number of solutions. 
Considering the restrictions for supersonic solutions, the demand can be made, that two possible solutions 
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(supersonic and the potential to shock down to subsonic flow) must exist. For two solution types it is necessary 
that two sign changes occur or that: 
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(2.83) 


The first condition is meaningful, the second condition is not permitted, since the minimum value for C 0 , 
occurs for C 0 > 2(y 2 -l)/y 2 . Hence, it is apparent that for two solutions to exist, which imply supersonic flow, 
with the potential to shock down to subsonic flow, we must bound Co, as follows: 2(y 2 - 1 )/y 2 <C 0 <2. The use of 
this relationship is a criterion for determining the feasibility of fully supersonic flow. It can be seen graphically 
that indeed this must be the case, by referring to Figure 2. 1 6. 

Fully, subsonic flow, on the other hand, is less encumbered. Here, it is necessary to demand only a single 
change in sign of the relevant polynomial. With this restriction it is possible to require that the first term: 
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C 0 >2 


(2.84) 


There is no upper bound on this branch, in that the subsonic branch is not well defined as M->0. 

Summarizing the results of this section, a family of restrictions for the local inversion polynomials, which 
are in turn related to the local flow field conditions have been presented, as presented in Table 2.3. 
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Flow Field Type 

(pu 2 +p) 2 

Restriction upon C 0 = — 

(pu)(puH) 

Mixed, supersonic/subsonic (i.e. transonic) 

2{y 2 -1) 

Comod^ 2 => subsonic 

7 

Fully supersonic 

2(y 2 -iy Y 2 <C 0 <2 

Fully subsonic 

C 0 ^2 


TABLE 2.3 Thermodynamic restrictions upon flow regime for the complete inversion relationship. 


2,4.4 Quasi- 1-d extensions 

Although the above analysis is exact for 1-d flow fields, it is necessary to correct the problem for quasi- 
one-d flows caused by variable area shrouds, such as shown in Figure 2.17: 
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h(x) 
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FIGURE 2.17 Variable area shroud geometry definitions. 


To perform this correction two strategies are employed: 

• Modification of the partial differential equation 

• local quasi- 1-d modifications to the local inversion scheme 

The first modification involves scaling the dependent variable <j> by h(x). Performing this substitution, the 
governing equation takes the form: 


dx dy h(x) 


(2.85) 


The new term added here has the same structure as a quasi- 1-d approximation. This term modifies the 
distribution of the conservative quantities, but it does nothing for local inversion process. Alternatively stated, 
the inversion from flux quantities to thermodynamic quantities is independent of changes of area. As an 
alternative to this situation, we consider the local control volume: 
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FIGURE 2.18 Variable area thermodynamic inversion method modification definitions. 


Our task is to define a local process which takes the unmodified primitive quantity (M, p, T, u, p) as it is 
predicted by the decoding algorithm (state 1), and then execute an isentropic variable area expansion or 
compression to yield (state 2). For example, consider the Mach number predicted to be M, using the decoding 
algorithm. Now, applying the Mach number area expansion using: 
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( 2 . 86 ) 


yields, M 2 the local Mach number with variable area shroud effects. Additional isentropic corrections are 
available for all static quantities using a similar control volume. Additional corrections are applied to the 
pressure term for very low Mach number flows, which are not uniformly retained by the isentropic model. This 
takes the form of modeling static pressure for low subsonic Mach numbers with a relationship of the form: 
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Which is a compressible Bernoulli’s type relationship. 
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2.5 Turbulence model 


The basic turbulence model developed for this research is an extended, algebraic model loosely based 
upon Prantl's second hypothesis for turbulent free jet and mixing layer type flows. There is stronger 
relationship with the model proposed by Reichardt, as described in Schlichting (1979). This is a form of the 
classical self similar analysis described by Tennekes and Lumley (1972). The range of extensions to this model 
is rather considerable. They include: 

( 1 ) Formulation in terms of flux quantities, rather than primitive variables 

(2) Near field analytical estimation of mixing constants and relationship to free shear layer mixing 
problems. 

(3) Analytical prediction of incompressible mixing layer growth rates using a non-linear stability 
analogy. 

(4) Semi-analytical estimation of compressibility effects upon mixing layer growth rates using a 
compressible vortex sheet model and non-linear stability analog, 

(5) Far field or boundary dominated mixing effects. 

(6) Vortical, 3-d, mixing enhancements due to the streamwise vorticity induced by lobed 
ejector/mixer nozzle configurations. 

As can be seen, the turbulence model for this flow is a considerable portion of the overall development 
problem. 

Tennekes and Lumley (1972) use the concept of self preservation or invariance to motivate their 
discussion on free shear layers (jets, mixing layers and wakes). The mathematical description for this type of 
“rule” would be termed “self similarity”. Our development proceeds from a somewhat different direction 
without explicitly invoking self-preservation or invariance. Nonetheless, the analysis developed in this study 
recovers these classical results, as indeed it must. When working with algebraic formulations, one of the few 
“fundamental principles” one can use to try to understand turbulent flow is self preservation or invariance. 

To avoid confusion, it should be emphasized that even though the turbulence model and governing mixing 
equations may be considered mathematically self-similar, the bounded internal flow that characterizes ejectors 
cannot be self similar due the existence of the boundaries themselves. This is the reason that it was not 
possible to solve the governing equations using some form of a similarity or Boltzmann (Bear(1972)) type 
variable. The fact that the local mixing behavior will be modeled using equations that are consistent with self 
similarity or self preservation is not at all inconsistent with this fact, since locally, the effect of the wall 
boundaries cannot be “sensed.” 

2.5.1 Two-dimensional mixing layer model 

The basic two-dimensional shear layer model is described in this section. As indicated the basic model is 
an algebraic formulation in terms of an effective viscosity. Starting with the basic "near field” (boundaries are 
not yet important in this region of the flow) differential equation, a relationship may be drawn between the 
effective viscosity and the thickness of the mixing layer. The result of this is the elementary turbulence model: 
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which is derived by analytically solving the partial differential equation: 
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for a mixing layer with a 0.99 mixing thickness definition, i.e. G(-8)-G|, G(8)-G 2 . It is possible to formalize 
the method by solving this equation: 
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Evaluating this equation at y=±5 and subtracting yields: 

0.99 (j) w = 0 + ~ <f> 2 o)erf(ri «>) 

-0.99 <{> 20 =1>-^(<P W ~ <j> 2Q )erf{Tl\) (2.91) 
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from which one derives: 

t]\ = erf~' (0.99) * 1.82 (2.92) 

Finally, decomposing, t]*, yields: 



(2.93) 


which provides the desired relationship between the effective viscosity and 8/x, the mixing layer spreading rate: 


kxe G = 



(2.94) 


It is clear from equation (2.94) that prediction of the turbulent viscosity is directly related to estimation (or 
empirical measurement) of the shear layer thickness, 8. A turbulent Reynolds number may also be written for 
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this problem. As indicated previously, the turbulent viscosity may be used to model heat transfer (energy 
transport) and mass transport by definition of appropriate turbulent Prandtl and Schmidt numbers (White 
(1991) or Williams (1988)). 

Recall, that the basic equations were perturbed using e=(Ui(rU 2 oy(Uio+U 2 o). The above analysis for the 
effective viscosity was based upon e G . Though, we expect these small parameters to be of the same order of 
magnitude, it is possible to write the effect upon the turbulence model, which accounts for the in the 
momentum equation (1 .22): 
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(2.96) 


It is not necessary to assume that the momentum scale and velocity scale ratio are 0(1), i.e. Zq!z » 0(1). Since 
for any problem, s G and e are both known, these scales are included in our analysis. 

Now for clarity, the closure approximations from Section 1 are repeated: 
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and the turbulent energy closure: 
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and the turbulent mass term: 
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Clearly, to complete the mixing model, an analysis is needed to predict (v efl /U)‘ and therefore 8/x. This 
development performed in the next section. 
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2.5. 1.1 Kelvin-Helmholtz instability 

Although the above equations relate spreading rate to the turbulent viscosity, it still relies on predictions of 
the shear layer spreading rate, 5/x. These predictions may be supplied using the classical, incompressible 
closure: 


, (UjzMA 
5 (u, + u 2 ) 


x 


(2.97) 


Physically, the spreading rate or mixing layer thickness may be defined by Figure 2.19. 



U 


FIGURE 2.19 Shear/mixing layer definitions. 


Though most work has developed this form using empirical or semi-empirical methods, here, a more 
formal theory based upon non-linear stability methods is developed here. We believe that this is viable, since 
the spreading rate for incompressible, turbulent mixing layers may be classically analyzed using dimensional 
reasoning and the concept of shear layer self-preservation, Tennekes and Lumley (1972), show that spreading 
rate, 8'=d8/dx, is a constant, hence, 8/x, is a constant as well. The goal of this section then is to predict, 8/x, for 
compressible and incompressible flows using a simple, yet analytically based method. 

It should be stated, however, that this work is not by any means the final say in predictive models for 2-d 
mixing layers. This field has very significant depth and breadth. Methods to study this problem range from 
basic dimensional reasoning to direct simulation. It is claimed, though, that this analysis provides a structured 
basis for predicting the form of the spreading rate, as well as, yielding deriving a compressibility correction to 
the above spreading rate equation. Additionally, the analytical basis of the prediction is desirable, in that, the 
ejector/mixer design problem is considerably more complex than the simple two-dimensional, incompressible, 
semi-infinite shear layer described by equation (2.94). As such, any method developed here must be 
generalizable to more complex configurations. 

A scaling model is developed for the growth rate of 2-D, compressible turbulent mixing layers. The basis 
of the development is the use of a time scale for growth of a dominant wavelength that is obtained using a 
nonlinear saturation time relationship with a temporal amplification factor obtained from quasi-unsteady, linear 
theory. The analysis extends to compressible flows by combining an unsteady, linear, vortex sheet 
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amplification factor and the same nonlinear theory. Results for compressible growth rates are in good 
agreement with experimental measurements and are expressed in terms of a convective Mach number. 

What is perhaps most striking about turbulent, free shear layers is the remarkable degree of repeatable, 
large scale structures exhibited in their growth and development. This structure is well documented in the work 
of Papamoschou and Roshko (1988) and Brown and Roshko (1974). Tam and Chen (1979); Morris, et al. 
(1990) and Liou and Morris (1992) have attributed the appearance of well-defined large scale or coherent 
structures in the fully developed turbulence to the relatively well understood, deterministic, growth of small 
initial disturbances, which are often modeled by linear and nonlinear small disturbance theory. They have 
connected instability wave theory to fully developed turbulence and their results are in reasonable agreement 
with the experimental data. However, these analyses are complex and require a numerical computation of 
linear disturbance growth rates from an Orr-Sommerfeld equation or the inviscid Rayleigh equation, as well as 
the computation of eigenfunctions or amplitudes from an auxiliary relationship, most usually some form of the 
turbulent kinetic energy equation. These wave solutions are then related to turbulence using a superposition 
technique or a hypothesized closure. The method used here relates to the previous works but with a simpler 
scaling relationship. It is shown that the result is a model for the turbulent spreading rate that agrees with the 
data and lends itself to a generalized approach. A similar direct relationship between small disturbance growth 
rates and turbulent spreading rates is discussed by Lu and Lele (1994), however, they only provide a limited 
justification and model development. 

The present model combines the temporal amplification factor (Drazin and Reid (1981)) and the concept 
of a nonlinear e-folding time (Youngs (1984), Andrews(1993)) to yield an estimate of the time required for the 
growth of the dominant wavelength. This temporal/spectral growth rate relationship is converted to a spatial 
estimate of the growth of mixing layer structures that determines the mixing layer spreading rate. The result of 
the analysis is a mathematical formulation for the mixing layer spreading rate that extends to compressible 
flows. 


2.5. 1.2 Linking linear stability with fully developed turbulence 

Introducing an analogy with Morris, et al. (1990) turbulent fluctuations are related to linear disturbance 
amplification factors, s and £ as: 
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The constant, p 0 , is to be determined later. The close relationship between the turbulent fluctuation closure in 
equation (2.97) and the linear disturbance theory developed in the next paragraph is deliberate. 

The amplification rate, s, for a 2-D, incompressible, vortex sheet mixing layer may be obtained using 
perturbation potentials (Drazin and Reed (1981)) as: 
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(2.99) 


(U i + Upki k(U 1 -U 2 ) 
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2 2 

Equation (2.99) shows that this flow is unconditionally unstable for any finite velocity difference. 

With the small disturbance amplification rate available, the sizes of the turbulent fluctuations are 
determined using Equation (2.96). The closure given by equation (2.97) has two unknowns: the amplitude 
constant p 0 and the wave number k. The wave number k has dimension 1/length, however, a free shear layer 
has no defined length dimension except for the mixing layer thickness, 8. Thus the following relationship for 
the wavenumber is proposed: 


2k 2k 
X ao<5 


( 2 . 100 ) 


where cto is a proportionality constant which is discussed later. 

Specification of the wavenumber using the shear layer thickness is a key assumption of the present 
analysis, and may be motivated by considering the approximate rapid deformation (Arpaci and Larsen (1984)) 
simplification of the turbulent kinetic energy equation: 


U—(-u'u')=-u'v'— ( 2 . 101 ) 

dx 2 8 


where AU=(U 1 -U 2 )/2. Since equation (2. 1 00) is in terms of the turbulent fluctuations, equations (2.97)-(2. 101) 
may be combined to yield an ordinary differential equation for 8(x): 


d_ 

dx 
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( 2 . 102 ) 


For analytical purposes, equation (2.102) is a complex nonlinear differential equation that cannot be 
directly solved. The complexity may be overcome if a suitable particular solution can be identified. From 
dimensional reasoning, the mixing layer width takes the functional form 8 oc x. If the functional form of 8 and 
k has been correctly postulated, equation (2.102) should collapse to an algebraic relationship between 
constants. Making the substitution: 


8=C S 



into equation (2.101), gives: 


Po ~ Cs 


(2.103) 


(2.104) 
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Thus, the turbulent fluctuation and wave number closures, equations (2.97) and (2.103), are particular solutions 
of the turbulent kinetic energy equation. Equation (2.103) indicates that the magnitude of the cross-stream 
fluctuation, p 0 , is related to the shear layer spreading rate. This is plausible because cross-stream fluctuations 
should scale according to the only available cross-stream length scale 5. 

The result in equation (2.104) indicates a strong connection between turbulent kinetic fluctuations and 
small disturbance stability theory. Further, the structure of the previous equations developed using dimensional 
considerations, has been shown to be correct. To close these equations, the constants Oo and C CT are evaluated in 
the next section. 


2.5. 1.3 Non-linear scaling model 

The previous section showed that the equations describing the shear layer spreading rate satisfy the 
turbulent kinetic energy equation if an appropriate wave number closure is applied. The dimensional reasoning 
that provided an appropriate closure for the wave number also suggests that the development of large scale, 
coherent structures in these flows does not require the solution of the previous equations. Instead, a model for 
shear layer growth rate is derived by recognizing that the fundamental time and length scales are contained 
within the linear theory and the overall convective frame of reference. The benefit of this type of analysis is 
that the spreading rate function can be obtained for a very general class of problems without the burden of 
directly solving a more complete set of equations, such as the boundary layer or Navier Stokes equations. 

Starting with the linear stability growth rate, equation (2.99), the e-folding time scale governing the growth 
of disturbances in the flow is l/s^ N multiples of the e- folding time, is the time required for the dominant 
wavelength, X d , to appear. A constant of N for different wavelengths implies a coherency between the 
developing large scale flow structures that arises through the saturation of the amplitude growth rate of 
instabilities. These concepts are related to the buoyancy driven mixing flow and jet break up models 
empirically developed by Youngs (1984) and Andrews (1993), respectively. From a survey of experimental 
data, they estimated that N* 9-10. 

The e-folding time is obtained from the real part of equation (2.99) as: 


L 2 

Sr kjdJrUi) 


(2.105) 


Similarities between supersonic and subsonic amplification factors mean it is not necessary to explicitly 
differentiate between supersonic and subsonic factors. Taking td=Ntd and further introducing the wave 
number definition k=2p/ld yields: 


L 


= t d U = 


al IN l 
2 n U1-U2 2 


(U, + Ui) 


(2.106) 


Introducing equation (2.105) into equation (2.106) and solving for 8/L yields: 

/ ^2S'=— 

L aN (Ui + U2) 


(2.107) 


Equation (2. 107) provides a formula for the spreading rate of a 2-D incompressible shear layer. 
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Although equation (2.107) represents a complete relationship for the spreading rate, the wave number 
closure constant, cto, has not been determined. No additional relationships are available to provide an exact 
determination of cto- An estimate of oto can be ascertained by considering the saturation length scale, L^L. 
This length scale was previously estimated as L=NU/s r to arrive at equation (2.106). If it is recognized that the 
available local length scale in a shear flow is 5, then equation (2.107) may be written: 

o \ J 

L W NS * — (Ui + Ui) (2.108) 

2k U1-U2 


Dividing both sides of equation (2. 1 08) by N5 yields: 


1 = 


go 1 

2 n £/i - Ui 


Oh + U2) 


(2.109) 


which should be considered an order of magnitude relationship rather than a strict equality. 
Now, if we recognize that the velocity ratio terms in equation (2.109) are bounded: 


JUi + U 2) „ 

< < oo 

Ui-Ui 


( 2 . 110 ) 


Thus, for equation (2. 1 09) to hold the term must be bounded by: 


0 < — <1 => 0<cc 0 <2k 

2k 


( 2 . 111 ) 


This is still a wide bound. However, oto is expected to be a universal constant, virtually independent of initial 
conditions and flow conditions, and so may also be estimated from experimental measurements. Additionally 
empirical measurements indicate that cto will be 0(1), so that the lower bound, i.e. cto = 0 will never occur. 
Following Andrews (1993), and inspection of the experimental photography of Brown and Roshko (1974), an 
average wavelength versus amplitude ratio of the largest scales was measured. This is done by measuring the 
streamwise length between any repeating portion of the large scale structures. This length is then scaled by the 
local mixing layer half width giving the ratio X/8. This ratio is observed to be approximately 4.0, hence from 
equation (5), cto=A/8 =4. This value is used in the current work and is within the upper and lower bound of 
equation (2.1 1 1). We note with some satisfaction that the arithmetic average between the upper and lower 
bound, oto= 7 t is within about 20% of the experimentally measured value oto =4. 

Equation (2.107) may be re-written: 


28 


lft (UrUj 

lN ±a ,,+u* 


( 2 . 112 ) 
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so that Papamoschou and Roshko’s (1988) constant density form of the convective velocity U=U c =(U 1 +U 2 )/2 is 
apparent. Equation (2.112) reproduces the functional form of the incompressible spreading rate functions 
given by Papamoschou and Roshko (1988), Brown and Roshko (1974) and Abramovich (1963) as: 


28 


U1-U2 
^(Ui + Uj 


(2.113) 


The result of equation (31), yields a prediction of C s =tc/ 2N. Table 2.4 favorably compares this analytical 
predicted value with other values. 


2.5. 1.4 Compressibility effects upon turbulent mixing 

In this section, the effect of compressibility on the spreading rate of the shear layer is discussed. It is 
our desire to use the same type of analytically based formulation as developed previously. Then an 
approximate extension the mixing layer relationship, which was used to relate the turbulent viscosity to the 
spreading rate, to handle turbulent flow. 


■ 

Analytical 

C 8 =7t/2N 

a 0 =4 

Brown and 
Roshko (1974) 

Papamoschou 

and 

Roshko (1988) 

Lange ly 
Curve 

C 8 N=10 

0.16 

0.19 

0.17 

0.14 

relative error 


17.3% 

7.6% 

14.3% 

C 5 ,N=9 

0.17 

0.19 

0.17 

0.14 

relative error 


10.5% 

0.0% 

21.4% 


TABLE 2.4 Proportionality constant, C 6 , for incompressible mixing layer spreading rate. 


To relax the limitations of the quasi-steady theory used to compute the mixing layer relationship, the 
complete unsteady potential equations are considered (Pai (1959)): 


M? f d 2 <fi, 

Ur dt 3 


+2u '^ ] =-ti 




(2.114) 


The form of the incompressible solution suggests a normal mode solution and interface disturbance of the form 
(Drazin and Reid (1981)): 
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<P, = <P J (y)e' k(x - u,)+sl 


(2.115) 




tk(x-Ut) MY 


A Galilean transformation has been made using the average velocity U=(U|+U 2 )/2. This transformation affects 
the phase speed, but not the temporal growth rate s. Following the method of normal modes the functions 
described by equation (2.1 15) are substituted into equation (2.1 14) and the resulting boundary value problem is 
solved. Applying far field boundary conditions, the interface kinematic relationship and the compressible form 
of the pressure matching condition yields the dispersion relationship: 


(s + ikAU) 2 (s-ikAU) 2 


aU T, 


a] <72 


; a, a k(l- — 2 (±Au-—f)i 

Gi K 


(2.116) 


where AU=(Ui -U 2 )/2. Miles (1958) solved equation (2.1 16) for “s” as: 


L > (U,-U 2 ) 

s = -aki[ Ml + 1 -(4 M 2 C + 1 )* T 2 , M c = — 


(2.117) 


We note, that this is a constant temperature solution. We may approximately relax the constant temperature 
restriction by introducing the correction term into (2.1 18) to yield: 


s = 



ki[M 2 c + 1-(4 M 2 c + l) 2 f 2 


Me = 


(UrUA 

2a 


(2.118) 


The explicit relation for the temporal growth rate, equation (2.1 17), is used to find the fully turbulent spreading 
rate using the same coherency concepts developed previously, giving: 


£_ n 

L ~ ° 2a 0 NU 


f(4 M 2 c + 1)2- Ml-1]2 


(2.119) 


Equation (2.1 19) is normalized by the incompressible solution, equation (2.1 13), to give a compressibility 
correction: 


S' 


M c 


( 2 . 120 ) 


This compressibility correction is for the spreading rate of 2-D mixing layers. It is dependent on the 
convective Mach number, M c . Later, it is shown to be consistent with the experimental work of Papamoschou 
and Roshko (1988) and others. 
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0 Papamoschou and Roshko ( 1 988), experimental data 
O BogdanofT ( 1 983), experimental data 

Unsteady solution, equation (2. 1 20) 

Generalized model solution; linearized growth rate from Gropengiesser (1970) 



FIGURE 2.20 Normalized spreading rate comparison of theory to experiments. 


Though, the formulation presented in Figure 2.20 is not valid for convective Mach number greater than V2 
a simple piecewise construction has proven to be effective. This has been the methodology chosen by us to 
correct for the effects of compressible flow on turbulent mixing. To extend the model for these convective 
Mach numbers, we employ a piecewise constant extension to apply the solution for Mc>V2. 

This analysis has done an adequate job in estimating the effect of compressibility upon the mixing layer 
thickness. This is the most important effect. A secondary effect of compressibility is to modify the mixing 
relationship used to relate the velocity mixing layer thickness to the momentum mixing layer thickness. 

It is appropriate to move away from the incompressible equation (2.94) and consider the more general 
equation: 
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d_ 

dy 
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d<j>^ 
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_ ^10 - <*20 
^10 + $ 20 


= + 


( 2 . 121 ) 
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In this relationship, it is necessary to include the effect of the variable density term on the left hand side. To do 
this a modification of the Lee -Dorodnitsyn transformation: 


\pdy 

o 

p 0 x 


( 2 . 122 ) 


is introduced, where y* is our “stretched variable y*=ey. Employing this new similarity variable, will modify 
the governing equation to take the form: 


d f p d<jf\ VA -i d<j>' 

-7-7- — -rr +77 w -7T 
drj \p o drj ) dr/ 


(2.123) 


Unfortunately, at this point it is not possible proceed without information concerning p/p 0 . This is formally 
possible to model, by solving the energy equation or applying some kind of particular energy integral (Crocco- 
Busemann). Instead, we will “freeze” the parameter p/p 0 . This term is then carried through the analysis and 
evaluated it at the appropriate time. Though this is by no means rigorous, it should capture the gross effects 
pressure gradient (at least better than linearizing at this point). Following this strategy, the above equation is 
solved to yield: 


1 + e erf 

[M 
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2 
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\Po J 



(2.124) 


Following the same strategy used in the incompressible development, the solution is matched to the 0.99% free 
stream (unmixed) conditions. This defines 5 the mixing layer half thickness. Applying this analysis yield the 
dimensionless number: 


77 *. = 1.82 



(2.125) 


Which clearly collapses back to the incompressible value for constant density flows. This is consistent with 
our previous formulation. 


2.5.5 Far field "boundary dominated" effects 

The previously described turbulence model is based upon free shear (semi-infinite) flow domains. Clearly, 
this type of model is valid for an internal flow for short shroud mixing problems. To extend the previously 


89 



developed concepts, in a rigorous manner the spreading rate function is modified by invoking the turbulence 
kinetic energy equation: 


<3c 

u ir v 


eff 


du 


3 

k 2 


— I -c„ — 


\dy 


l 


(2.126) 


where the terms refer to: [transport]=[production]-[dissipation] (Arpaci and Larsen (1984)). By introducing the 
approximate closures: k=(Au8*) 2 , 8*=8/H , v turt) = 8Au and 1=H, we may write: 


dS 2 

— - Cs e(l - S 2 ) 


<?_ tanh (j£gx) 
x x 


(2.127) 


The solution of this elementary Riccati equation (Boyce and DiPrima (1965)) provides a limiting value to 
growth of the mixing layer thickness, i.e. 5 max <l . 

One limitation to this formulation, however, is that as x-*oo, the 8/x->0. This is an artifact of the 
derivation, since v cfT /U 0 was derived based upon a near field (no wall influence relationship), i.e. equation 
(2.89). As such, we cannot expect the turbulence model developed here to be valid for very large values. It is 
possible to estimate the estimate the size of these large values by considering the streamwise value for which 
tanh(C 5 x max )=0.99. Solving for x max yields: 


x 


max 


= — tanh-’(0.99)»l 
C$s 


(2.128) 


Since this study is concerned principally with short shroud aerodynamic ejectors and the acceptable streamwise 
length is on the order of 15/e, this is seen to be a minor limitation. 


2.5.6 Vortical mixing enhancements 

Mixing within the ejector shroud is the key to satisfying noise constraints for high speed civil transport 
problems. Obtaining this mixing within a sufficiently short streamwise distance, which implies a short and 
therefore small weight penalty, is a key technology for the high speed civil transport program. Methods to 
obtain this goal involve the distortion of the splitter plate into a form that enhances mixing. This type of 
corrugated splitter plate is illustrated in the Figure 2.21 : 
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FIGURE 2.21 General vortical flow field definitions. 


For a practical vortical mixer, these corrugations are greatly exaggerated and are often referred to as chutes. 

The basic mixing enhancement mechanism caused by the introduction of convolutions is the formation of 
a rotational or vortical component. The axis of the resultant velocity vector is in the streamwise direction (or 
“x” direction); hence the rotational distortion is referred to as streamwise vorticity. The effect upon the fluid of 
the streamwise vorticity is the dramatic elongation of the interface i.e. mixing layer, which separates the two 
flows. This elongation and ultimate distortion of the mixing layer interface dramatically increases the 
“effective” two-dimensional mixing layer. 

Since the model developed here is a strictly two dimensional approximation the introduction of this 
inclusion of the inherently three dimensional effects of vortical mixing would seem to be beyond the reach of 
this model. Fortunately, though, the three dimensional effects of vortical mixing can be mapped into an 
effective 2-d mixing layer. Our strategy to do this, is to model and predict the three dimensional interface 
elongation and distortion, then to use this solution to enhance the 2-d shear layer. Effectively, we are 
decoupling the 2-d mixing layer from the three-dimensional roll up. Though only a hypothesis, it seems to be a 
valid approach as discussed by Qui (1996), Fung (1996) and Marble (1985). Additionally, to introduce a 
model of greater complexity, would be to invalidate the scope of our current work which is a preliminary 
design tool. 

The task, then is to develop a relatively simple model which predicts the behavior of the kinematic 
interface (or slipline) between the two mixing stream. Clearly, this roll up or twisting cannot proceed 
indefinitely. Ultimately, the interface must interfere with itself and with adjacent (and counter-rotating vortical 
cells. At this point, turbulent diffusive effects dominate and the roll up structures collapse onto themselves. 
This increase of interface, followed by a collapse mechanism, brings to mind an analogy with water waves 
breaking on a beach. The water wave surface or interface grows in an inviscid manner. Eventually, though, 
this growth steepens to a point where the wave can no longer hold itself, (and actually would become double 
valued). To prevent this double valued effect it is necessary introduce the concept of a normal shock or 
discontinuity. In reality of course, this breaking wave and discontinuity cannot occur without local smoothing 
due to viscosity. The mechanism described is one classically modeled by a Burgers equation (Anderson et al. 
(1984) or Logan (1987)). 
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Our modeling of the vortical mixing phenomenon uses the above breaking wave solution as an analogy. 
Proceeding, the governing equation for the interface problem is written (Drazin and Reid (1981)): 


cos^=t/,^-sin^^ 
ax az 


(2.129) 


where V 9 is the local rotational velocity caused by the streamwise vorticity. Geometry for this flow is shown in 
Figure 2.22. 



FIGURE 2.22 Vortical flow interface distortion definitions. 

Now, very early in the flow, i.e. close to the splitter plate, where 0«1, the kinematic equation reduces to: 


v.u£ 

y ff U 0 -a 

ax 


(2.130) 


Integrating and replacing V e by an inviscid, point vortex solution, V 0 =r/27tr, yields: 
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near 


i r 

U 0 2nr X 


(2.131) 


which is valid in the “near field”. This relationship is approximate, in that near field roll up does not act like a 
simple vortex, hence we generalize: 
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inv^vort 
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2nr 


(2.132) 


where C^von is a semi-empirical constant which takes the approximate value of about two. 

Farther downstream, where interference effects will begin to dominate, it is no longer possible to utilize 
the near field solution. Here, we make a far field approximation valid for n/2<6<2n: 
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(2.133) 


Substituting into the interface equation yields: 


■—V =U — + — Vg — 

v 6,i u e * , T v 0,t o 

z n dx z n dz 


(2.134) 


This equation is a damped non-linear wave (Logan (1978)). It is necessary to be able to model the left hand 
side of this equation. To do this a Reynolds stress analogy is applied, yielding: 
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effort 
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dz 1 


(2.135) 


The effective, vortical kinematic viscosity will be estimated later. Accepting this new closure, equation (2.129) 
may be rewritten: 


U K + ±v 
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eff ,vort 
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(2.136) 


This is a classical Burgers equation. 

Solution of the Burgers equation, given sufficiently simple boundary conditions turns out to be 
straightforward. To illustrate this, we consider a canonical problem of the form: 


du du 
— + u — = 


dt dx 



(2.137) 


where, u, represents any general scalar variable. Now if we have the boundary conditions: 
u(-oo)=u war and u(oo)=Uf ar „ presented in Figure 2.23. 
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FIGURE 2.23 Breaking wave analogy model for vortical flow definitions. 


It may be noted, that for these boundary conditions to make sense for our physical problem, that we must 
translate the problem downstream some amount x c . See Figure 2.23. With these conditions, a traveling wave 
solution is possible, hence we search for a traveling wave similarity solution of the form: 


u(x,t) = f{t-kx) = f(s ) 


(2.138) 


where k is computed as a part of the solution. Substituting yields the ordinary differential equation: 


f-kfT= v* 2 /" 


(2.139) 


Integrating this equation twice and applying the associated boundaiy conditions yields: 


u = 


^ near ^ far** 
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(“near ~ ^ far) 


2 vk 


(2.140) 


From the far field conditions, we obtain, k=2/(u near +u far ). 

From the above equation, we can return to the actual physical problem by noting, 
x-x^Zo/Vej z, t=x/U 0 , v=v cff)VOf1 (z</Ve ,) 2 and u=£. With these conditions the above equation becomes: 
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V effort* 2 * 
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(2.141) 


It is helpful to associate the unknown Xc, with a “shock location,” a location of rapid collapse of our vortical 
system. The shock location is found by noting that at x=0, £=0.99£ near . It is worth noting the similarity with a 
99% boundary layer thickness definition. If this is the case, it is possible to write the non-linear equation: 
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Which is a non-linear equation for Xc. 

To proceed, it is necessary to associate the near field £, near constant, with our near field solution: 


£ 


near 


i r 

U Q 2 nr 


(2.143) 


It is worth commenting that although the above analysis required a constant near field solution, we have 
approximately extended the analysis to include the variable function. If the previous analysis is to be suitable, 
it is certainly not possible to merely substitute the previous equation into our shock layer computations. Instead 
it is necessary to define a maximum value for equation (2.413). This maximum is related to the early linear 
growth phase, followed by the shock decay. This maximum is denoted as £ max . 
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(2.144) 


It is necessary to further, modify the initial interface function to include the additional surface area 
associated with the corrugations of the mixing layer. This will occur even in the absence of vortical mixing, i.e. 
r=0. This initial interface length will ultimately be damped out by viscous/turbulent mixing. To account for 
this effect we include an exponential decay term, which scales according to turbulent mixing only. Writing this 
additional interface length term: 


1 Vtyrbf 

( O )*' 1 14 ® < 2 - l45 > 

Finally, the far field solution, is characterized by no stretching of the interface, hence ^ f „=0. Clearly, this will 
simplify the previous equations dramatically. 

The above relationships contain the unknown £, max . This unknown is associated with our critical location 
formula, equation (2.144). It is now necessary to compute £ max , and x c . This may be achieved if we consider 
the shock location and breaking phenomenon, to be a patching between an initial solution (our near field 
equation) and the far field. For a steep shock front, this is a reasonable approximation. Mathematically, this 
patching takes the form: 
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max 


= C„ 



(2.146) 


where equation (2.144) has been rewritten: 
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(2.147) 
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Solving for £ max yields: 


= 18.38C, 


V effyor, Z <> 


inv t vort 


v 2 e,i 


(2.148) 


To complete the solution, it is necessary to model the effective vortical viscosity term. We choose to use 
the average velocity and the initial lobe heights as relevant length and velocity scales; v eJtv0(tiail =CoUoZo. C 0 is a 
semi-empirical constant related to the number of vortical rotations that we might expect before collapse of the 
defined interface structure. It’s value is approximately 16.-32. With this closure, it is possible to write the 
maximum interface height in a standard non-dimensional form: 
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(2.149) 


This non-dimensional form is significant, in that this relationship between interface maximum distortion and 
vortical Reynolds number (Rev^r/Uo^), has been observed from the full Navier-Stokes, computational 
experiments of Elliott et al. (1992) and Fung (1995), as the appropriate scaling law for strong vorticity flows. 
We have grossly reproduced a composite of their figures, to exemplify this relationship. This is presented in 
Figure 2.24: 


log 10 (mixing 
augmentation rate) 



log 10 (Re vort ) 


FIGURE 2.24 Scaling law from the Navier Stokes computations of Qui as reported by Elliott et al. (1992) 

and Fung (1995). 


In our work, of course, the mixing augmentation rate is associated with Hence from equation 

(2.149) it is apparent that we have reproduced the high vortical Reynolds number regime functional 
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relationship, i.e. the power of -1/3. It is gratifying, that our coarse, analogy based model has recovered this 
dependence. 


At this point, our development of the function describing the shape of the interface separating the two 
streams is complete. It is necessary, to estimate the surface area of an elemental strip of the vortically mixed 
interface surface. This is formally done by computing an integral of the form: 


f-Ji + [(«*•*»’ ] i<fe 


(2.150) 


Since the interface function is essential only a function of x, the integral may be approximated on a per lobe 
basis: 
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(2.151) 


As a summary, it is desirable to collect the pieces of our analysis which consist of the. 

• -the interface growth and decay 

• -shock location and initial maximum 

• -and the initial interface function. 

Combining one obtains the final relationship: 
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(2.152) 


where: 

£o(x) = Additional surface area length due to splitter plate corrugations 

X - lobe wavelength 
r = streamwise vorticity circulation 
ho = total height of offset 
r =70 

'^vort.iriv 

as, shown in Figure 2.25. The initial wave length relationship takes the form: 
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(2.153) 
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(2.154) 


which may be recognized to be an incompressible simplification of the mixing layer growth rate equations. 



FIGURE 2.25 Vortical flow geometry input parameter definitions. 


Finally, to close the above relationships, the vortical circulation, T, or more importantly it’s non-dimensional 
form, must be estimated. The vortical Reynolds number Revo,, is related to the chute geometry using the 
inviscid/continuity relationship of Skebe and Barber (1988): 
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(2.155) 


This is the major effect of vortical mixing. We also hypothesize an additional effect upon the 2-d shear layer 
thickness 8. It is expected that the addition of streamwise vorticity in addition to causing greater surface area 
will be to increase the spanwise vortical effects. The “equivalence” of this effect to the streamwise vorticity is 
motivated by considering a 3-d vortex sheet equation: 
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(2.156) 
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To compute the unsteady growth rate as in Section 2.5.4 it is necessary to introduce a normal mode solution 
that takes the form: 


<f>. = (y) e ik < *+k-uo*st g = e ik(x + iz-uo+s, 


( 2 . 157 ) 


Where 1, is a dimensionless wave number for disturbances in the z direction. At this point, one can perform 
virtually the same analysis as in Section 2.5.4 to obtain: 

[s + /ft(A£/+AK/)] 2 _ \s-ik(MJ+AV t l)] 2 (2158) 

a] a i a] a 2 

From this equation, the “equivalence” of disturbances due to shear, whether due to the mixing layer, AU, or due 
to vortical effects, AV Z is apparent at least within the scope of the vortex sheet formulation. 

Physically, it is possible to use a schematic of a single vortex and the hypothesized equivalent forms to 
help motivate the previous mathematical structure as shown in Figure 2.26: 



Using the development, we estimate V z =2U ave Re vorl . This approximation is used to both modify the convective 
Mach number, as well as modify the expected closure for 5. Finally, we state that this modification goes along 
with at least the authors intuition, that vortical enhancement is not merely the wrapping and distorting of a 2-d 
mixing layer sheet, but that the growth rate of the sheet, i.e. 8(x,z), is modified by the vortical enhancement. 


2 . 6 Wall friction and heat transfer 

To develop a complete simulation tool the effects of friction and heat transfer are included in the analysis 
in an approximate manner. This takes the form of assessing compressible, turbulent skin friction and (through 
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Reynolds analogy) heat transfer losses. These are reflected in the governing equations through a wall loss 
term: 


gw _ +p) wall _ 2 

dy dy h(x) 


(2.159) 


The skin friction closure is compute through a locally fully-developed, locally adiabatic, internal skin friction 
relationship developed by De Chant (1993). This implicit relationship takes the form: 


V2 arcsin(l-a 2 ) ,/2 1 . 1 (^2 


fiMl 


a 


= -ln(Re,C“) + - 

AT v } K 


. _ , -( 1 + 2 ®)^ 
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(2.160) 
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l + 0.3(-)Re A C; 2 (l-a 2 )^ 


where: 


a 


2 




(2.161) 


and Re h -Uh/v wa n, B-5.5, k-0.4, cd-0.76, the power law exponent for the air viscosity and r=recovery factor 
and for turbulent flow, this takes the form: r=Pr 5/3 (White (1974)). 

2. 7 Theoretical conclusions 

This section has described and developed the mathematical models used to provide an ejector mixer- 
nozzle analytical tool. At this point, it is worth summarizing. Section 1, provided the basic conservation 
equations using formal perturbations methods. Section 2, developed the initial condition models, a combined 
numerical/analytical decomposition solution method and the necessary turbulent closure. Coming up are 
Section 3, discussing numerical methods, and Section 4, which outlines “in house” experiments associated with 
the ejector. Finally, Section 5 compares the model to experimental results and provides conclusions and 
recommendations. 
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3. NUMERICAL METHODS 


Section 3 discusses the numerical methods used to solve the governing equations developed in the 
previous section. Since a heavy emphasis has been placed upon analytical developments, these models are 
relatively simple. In spite of this simplicity there are some relevant and interesting here, especially, the higher 
order differential equations solver. 


3. 1 Introduction and motivation 

Implicit finite difference relationships and their use in solving 1-d parabolic partial differential equations 
are developed and discussed. Differencing methods of this form have high accuracy in terms of truncation 
error, while requiring limited support. This type of formulation has been applied to a high efficiency, 
combined analytical/numerical fluid flow model. Since the flow model has been derived such that a relatively 
simple, linear, scalar parabolic equation governs the flow; it is appropriate to develop the numerical methods 
for a general, scalar heat equation. Methods to resolve initial condition singularities in the physical problem 
are discussed. The streamwise marching portion of this problem is differenced using both Crank-Nicolson and 
a three point backward fully implicit method. Stability of the implicit finite difference/Crank-Nicolson method 
is analyzed using von Nuemann’s method. Heuristic extensions are drawn for the three point backward 
approach. A brief discussion of the method for inversion of the linear system characterizing the implicit 
difference methods is also included. Finally a comparison of strengths and weaknesses and recommendations 
for further research for implicit differencing method are presented. 

Hirsch (1988) defines implicit finite difference formulas as expressions where derivatives at different mesh 
points appear simultaneously. Relationships of this form have high accuracy in terms of truncation error, while 
only requiring limited support. For this reason, this type of formulation is rather attractive for our actual 
research problem, namely, a high efficiency, combined analytical/numerical fluid flow model. The model has 
been derived such that a relatively simple, linear, scalar heat equation governs the flow. It is then appropriate 
to develop the numerical methods for a scalar heat equation. The marching portion of this problem is 
differenced using relatively standard techniques. Both Crank-Nicolson and a three point backward fully 
implicit method have been employed. These methods have a lower order truncation error, but are A-stable as 
defined by Ferziger (1981). A further rational beyond stability concerns for the combined approach of high 
accuracy differencing in the cross stream (spatial) coordinate and lower accuracy in the streamwise (effective 
temporal) coordinate is that the physical problem has very rapid changes in the cross-stream direction (the 
actual initial condition is modeled as a Heaviside step function) and relatively gentle changes in the streamwise 
direction. 

An outline of our discussion: 

• Derivation and general discussion of a variable step, compact, implicit finite difference algorithm applied 
to a general, linear heat equation. Relationships to other methods are also considered. 

• Resolution of singular behavior: tension spline bases and analytical decomposition. 

• Conservative (finite volume) implicit differences. 

• Streamwise differencing and von Neumann stability analysis of a heat equation. 

• Block tridiagonal matrix solution algorithm. 

Having described the numerical methodology, some recommendations on current and future uses of these 
techniques as wazzu compared to other methods are made. 

3.2 A variable step compact implicit finite difference algorithm 

As indicated in the introduction, Hirsch (1988) defines implicit finite difference formulas as expressions where 
derivatives at different mesh points appear simultaneously. These methods have high accuracy with limited, 
here three point, support. Many methods and associated terminologies have be used to derive implicit formulas 
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for derivatives, including compact methods Kreiss (1973) and Hirsch (1975), spline methods, Rubin and 
Khosla (1977), and Mehrstellen or Hermitian for their relationship to Hermitian finite elements, Collatz (1966). 
As pointed out by Hirsch (1988) or Peyret and Taylor (1984), formulas of this type may be systematically 
developed using Taylor series methods. This approach is followed, except it is generalized to variable grid 
spacing. 

Considering only the right hand side portion of the linear, 1-d, constant coefficient, parabolic partial differential 
equation 


cfy_ <?$ 
dc a dy 1 


d<b 

+ b— + cS + d 

% 


(3.1) 


implicit formulas for the “y” derivative terms on a variable grid spacing are developed for some function f 
(where clearly Referring to Figure 3.1 the variable grid is defined. 



f. 

J 



j-1 


j 


j+1 


< 

A. 

1 


>< 



> 


FIGURE 3.1 Variable grid definitions. 


Beginning by writing the Taylor series expansion for f 


-^A V < V^aV“ > , + 0(A)’ 


(3.2) 


/ "y+i = f"j+Aj*if"'j+Tj;A 2 j+if W j + ^A 3 /+l / (5) , +^-AVi/ (6) y + 0(A) 5 (3.3) 


Combining terms in equations (3.2) and (3.3), yields: 
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(3.4) 


= /' /+ Ia , A, +iA J A / „(A / „ - A + 


A,.,+A, 


+0(A)* 

Now, it is desired to eliminate the second order term, to yield a higher order relationship. Consider the 
function f, expanded in Taylor series: 


/>-, = f, - V,4aVV~ a V + 


(3.5) 


+0(A) 6 


- f, + V,4aV, 4< 4’-r^4V‘ > / + ^A’,/”, + 


(3.6) 


+ 0 ( A) 6 


Equations (3.5) and (3.6) are combined to yield the term: 


A’y/,,1 - ( A’/ + AVi)A/, + A 1 ,/,-! 2 (a 2 ,-AVi) ._L„ „ , TO . 

A : ;+1 A 2 j(A j+1 + A ; ) ■ A A „ _/+/ ' + 

(3.7) 

+iA / A / „(A„,-A,)/« ! », + 0(A)' 


and rewriting equation (3.4) to eliminate the second order term in (3.7): 


_ a,)/’>, + 


+0(A) 4 


A ; + i + A , 


(3-8) 


combining terms and simplifying yields the implicit relationship for the second derivatives: 
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(3.9) 


^h.,r,- l +5(A,., + A,)r J +A J r,. 1 ] 

A’;/;., - (a 1 , + AVl]/; + A 1 ;/;., | 1 
A 2 7 +iA 2 7 180 


, (a 2 ; — A 2 j+|) 

A A>. 

AyA y+I (A y+1 - Aj)f^ 5) j + 0(A) 4 


Equation (3.9) is our implicit relationship for the second derivatives for a variable grid. It involves function 
values as well as first and second derivatives. This structure poses no major problem, since the solution 
methods presented here do not require an explicit formulation, but will solve the resulting implicit system 
simultaneously. The accuracy of equation (3.9) is formally third order for variable grid problems, but for 
gentle grid changes acts more like a fourth order method. This same claim is made for variable space problems 
(though explicit) by Ferziger (1981) and is consistent with our experience. For constant grid spacings, the 
classical implicit formula for the second derivatives is obtained: 

^■[/"y-.+10/'>/" 7+ i] = — 2 A { y+/> "' +0(A) 4 (3.10) 


In an analogous manner, the variable grid spacing implicit formula for the first derivatives is derived, to yield: 


^[A y .,/y,+2(A y . 1+ A,)/ J+ A J /’,, l ] 

+0(A)' 


a^i^-aV^-aV;., 

2A , +1 A, 


(3.11) 


which also reduces to the expected constant grid first derivative implicit formula: 


= / y > 1 2 A /; " 1 +0(A)4 < 3I2 > 

Equations (3.9) and (3.1 1) are the implicit relationships used in our analysis. 

As indicated, they are a of an overall system. This system may be more clearly discussed by introducing 
the new variables: <|> functional value, F=first derivative and S=second derivative. Hence, re-writing equations 
(3.1), (3.9) and (3.1 1) in terms of these new variable yields the system: 


dx. 


- aSj + bFj + c(j>j + d 


(3.13) 
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(3.14) 


1 (A 2 ; - A 2 7 + l) 

+ 5CA,., + A ft + V ; „] 4- — > ■ F, = 

- A ^-' - (A ’ + 4 a A .,<a,„ - A ^’’ * 0(A)< 

A A j 1 oU 


r 1 >$}+ 1 + j ~ ^ J +l Vy - ^J0j - 1 

[A, +1 ^-, +2(A >+1 + A y )/5 + A^,] = ^ 2A ^ |A ; i " 


(3.15) 


+0(A) 4 


It may be noted that the streamwise derivative, i.e. o/dx has not been treated yet. A discussion of this 
discretization is postponed until later. 

From equations (3.13), (3.14) and (3.15), the three equations needed to close the three unknown functions 
are available. Further, at any marching (streamwise) plane the matrix structure of the system is block 
tridiagonal since the support was three point. Inversion of block tridiagonal systems is described in Section 
3.5. 


The above governing difference equations are applied to any interior point within the flow. Appropriate 
boundary conditions are required. Following Hirsch (1975) the question is posed: what is the maximum 
accuracy one can ask for at any boundary point (necessarily involving only two adjacent grid points) and any 
combination of function, first and second derivatives? The answer is the relationship: 



with an analogous relationship at the other boundary. Equation (3.16) is also obtained using splines by Rubin 
and Khosla (1977). Equation (3.16) supplemented by equation (3.13) and the boundary type: e.g. 4>j = 0 for 
Dirichlet or Fj=0 for Nuemann boundary conditions. (These equations are the required three equations in three 
unknowns at the wall node). The benefit of this closure is that it avoids the boundary condition problem faced 
by five point support difference techniques, which typically must resort to a locally lower order method near 
the boundary. 

Realizing, that partial differential equations may be analyzed using many different methods, it would be 
useful to summarize alternative techniques. These are presented in Table 3.1. The basis function is used as a 
key parameter separating them. Comments on the strengths and weaknesses of methods are provided in the 
conclusions: 
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Method 

Basis Function 

finite difference 

polynomial (local) 

implicit finite difference 

polynomial in terms of function and derivatives; (global) 

traditional Galerkin 

polynomial (global) 

finite element 

polynomial (local) 

spectral 

trigonometric or Chebyshev (global) 


TABLE 3.1 Numerical solution methods for partial differential equations. 


3.3 Resolution of singular behavior : tension spline bases and analytical decomposition 

In the following sections, the initial condition singularity that was developed previously is now discussed 
with respect to differencing methods. Numerical methods, which are intended to better resolve the rapid 
changes in the initial condition, are discussed. Then the numerical/analytical decomposition is described for 
completeness. 


3.3.1 Flow singularities and analytical methods 


The idea behind deriving implicit equations based upon tension splines was motivated by the singular 
nature of the initial condition for the physical problem that is fundamental to this project. The step function 
initial condition is shown in Figure 3.2: 


wall shroud 



centerline 


FIGURE 3.2 Schematic of flow initial conditions. 


Additionally, governing equations for this problem may be approximately written for the conservation 
quantities in terms of the general linear parabolic equation: 



&± 

dy 2 


(3.17) 
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with the boundary conditions: 


with the initial condition: 


and $ is defined by: 


0) _ d<fi(x, l ) 
dy dy 


</>(0,y) = 
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< >- < 

h s <y<lj 
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<H x >y) = 
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pu‘ + P 

puH 

pn - 


(3.18) 


(3.19) 


(3.20) 


which are the flow quantities. These are, of course, the canonical forms of our mixing equations described in 
section 2. Here they may be considered any scalar simple equation. As before, it will be appropriate to make 
the local assumption, that the function a(x) may be approximated by the linear function a(x)=a x. 

As illustrated the flow is discontinuous at the interface between the primary and secondary streams. It was 
felt that a tension spline might provide a more satisfactory resolution of the initial condition, since tension 
splines are purposely designed to zones of large curvature. In fact for a very large tensioning parameter, a, the 
tension spline acts like a linear (i.e. “connect the dots 11 ) basis. Unfortunately, though this spline showed 
improvement, it still did not resolve the input condition adequately. Further, large values for a, cause 
degradation of the truncation error order. 

As indicated in Section 2 the poor performance of the strictly numerical model may be analyzed further. 
Upon considering the implicit relationship governing the second derivatives, equation (3.9), the error term is 
found to be: 


error = " A)/'"' + OW‘ (321) 

which is formally fourth order for constant grid spacing. Unfortunately for very rapid changes in f, the high 
order derivatives, such as, f< 5) , i* 6) are very large. In fact they may be as large as 0(1/A) 4 which would yield a 
very poor solution with errors 0(1). 

Further, this deficiency was shown clearly when the divergence theorem constraint, i.e. mass conservation 
which demands: 


l' o 0(x,y)dy = 1 


(3.22) 
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was considered. Unfortunately, this constraint was not properly honored by the strictly numerical solver due to 
the singular behavior near the step input. Physically, failure to observe the mass constraint is not acceptable. 
Ultimately, it became clear that the only possible way to appropriately deal with this singularity was introduce a 
local analytical solution that models the discontinuity. Examples of the use of basis or trial function (Galerkin 
terminology) come from Fletcher (1984) as well as others. Other examples of local treatment of rapidly 
varying functions near singularities include: steady state well functions used to resolve local radial well 
behavior in petroleum field simulation, Ewing (1983) and the use of “wall functions” to resolve local, near wall 
shear stress and velocity, for turbulent, viscous flow, Anderson et al. (1984). In the problem considered, here 
instead of developing a local special differencing method and blending it back into the overall system , the 
linearity of the governing equation itself was used to perform a global decomposition. 

This global decomposition begins by recognizing the requirement for an analysis capable of dealing with 
extremely rapidly changing functions. Fortunately a solution method from classical analysis is available which 
is based upon distribution theory rather than continuous functions, described previously in Section 2 and 
repeated for convenience. 


1 00 
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(3.23) 


It is important to note, that for x«l this relationship recovers the step input. In passing it may be noted, that a 
straight forward eigenfunction solution may be written: 


_ 2 00 i x : 
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(3.24) 


where: 


0=0 io h s + (l-h s )0 2o = l 


As pointed out previously, Equation (3.24) is formally an exact solution to our governing equation. 
Unfortunately, this equation is very slow to converge for x«l, 0(l/n), and suffers from Gibb’s jump 
departures. Worse, x«l is precisely where it is desired to obtain good converge since this is a region of great 
interest in the physical problem. 

Although equation (3.23) is exact, and does not suffer from the near field (small "x") limitations that the 
eigenfunction expansion solution, equation (3.23), did, it is still in the form of an infinite series. Note that in 
the near field, the solution converges extremely quickly. What is proposed and actually implemented is a 
combined numerical and analytical solution of the problem. Since it is linear this is easily effected. Consider 
the decomposition: 


0(x, y) = 0 an (x, y) + 0 num (x, y) 


(3.25) 


with the governing equation: 
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with the new boundary conditions: 


Wmm 

dx 


(3.26) 


/ i 9/1 um $ • dnum 

= a(x) . — ~ a x 


dy 2 


8y 2 


d n um 0) dan 0) dnum (•*•> 1) _ 9pn 1) 

dy dy dy dy 


(3.27) 


and the new initial condition: 


tum(O.y) = 0 


(3.28) 


The overall solution is then, the sum of the two functions. 

It may be noted, that a “patched” combination of the Green’s function expansion solution, equation (3.23) 
(useful in the near field), and the eigenfunction expansion solution (3.24) might provide a useful and efficient 
solution method valid in both near and far fields. This question is somewhat irrelevant for the 1-d scalar 
problems described here, since the numerical portion of numerical/analytical global decomposition is relatively 
inexpensive. On the other hand, when considering a 3-d mixing scheme (See Section 5.), numerical inversion 
may become significantly more expensive, thus making fully analytical solution more competitive. 


3.4 Streamwise differencing and von Neumann stability analysis 

The previous derivations have concentrated on cross-stream, “y” direction differencing. In this section the 
streamwise “x” discretization is discussed. Two methods, a completely implicit three point difference operator 
and a two point, semi-implicit, Crank-Nicolson operator, are considered. Then the more restrictive operator 
(the two point, Crank-Nicolson method) is tested for stability using the von Neumann method. 


3.4.1 Streamwise differencing 

Two different streamwise discretization methods are discussed. Consider first; the computational 
molecule for the second order, 3 point backward, fully implicit operator (which was ultimately the method of 
choice) presented in Figure 3.4: 
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FIGURE 3.3 Computational molecule for the three point backward discretization method. 

From Figure 3.4 it should be noted, that the “i” level is where the viscous, right-hand side of the governing 
equation, equation (3.1) is written. Figure 3.4 indicates that a variable grid, variable in the marching or “x” 
coordinate has been used. In a well-known way (we actually used Taylor series), the variable step formula for 
the stream wise derivative may be written: 
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2A,+A,., 
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(3.29) 


and for constant grid spacing recovers the 3-point backward formula: 

(3.30) 

With equation (3.29), the equation set (3.13), (3.14) and (3.15) is fully discretized. Rewriting (3.13) with 
double subscripts (i=streamwise, j=cross-stream). 


dx 
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(3.31) 


= aS iJ +bF iJ +c<f> u +d 


This is completes the computational set actually used in this work. One objection to the three point 
backward formulation is that it is not self starting, since 2 levels of information are required to begin the 
marching procedure. Actually, this is not a major difficulty, since the near field expansion equation (3.23) is 
very accurate in the near field even for a small number of terms. 

From the above derivations, it is apparent, that the method that is applied has a truncation error 0(A 2 * A 3 ’ 
4 j). Our motivation to use a lower order method is a consideration of stability limitations. A multistep, 
(Adams, Adams-Moulton) (Ferziger 1981) or predictor corrector method, i.e. of the Runge-Kutta methods 
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could have been used to gain higher marching accuracy. It was felt, though, that this trade of higher order 
accuracy would come at an unacceptable limitation in stability. 

3.5 Block tridiagonal matrix solution algorithm 


The implicit method that has been described requires the solution of a system of block (at minimum [3x3], 
and possibly more if there is a system of equations or non-linearity) tridiagonal linear equations. The matrix 
structure for this type of system may be conveniently described: 
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(3.32) 


where the [A], [B] and [C] terms are [3x3] submatrices and the x and d terms are [3x1] vectors. The x vector 
would be associated: x=[<j>,F,S] T . This matrix may be inverted in an efficient manner. Routines are provided 
by many workers, for example, Anderson et al. (1984). 

A simple way to understand the block tridiagonal algorithm is to merely generalize the scalar LU 
factorization method. This generalization is described by: 


[5] i =[5^-[^]'{[5r , }" , [c^ , 


(3.33) 


d‘ = d*-' - [A]{[Br i y : d‘-' (3.34) 

for i=l,2,3...N. At i-N the back substitution is begun: 

*-{[«]") V (3.35) 

and : 

x'={[B]fV-[C]V*') (3.36) 


where i=N-l ,N-2,... 1 . 


The above is a summary of the algorithm that are used to invert the block system described by equation 
(3.32). It should be noted that we would never compute the inverse of the partitioned matrices such as 
described in the above algorithm, instead one would always set up an equivalent problem, i.e.: 


Ill 



(3.37) 


* = {[*]"} V « [B) N x N = d N 


where the linear system would be solved using Gauss elimination with partial pivoting. 

This completes the summary of our block tridiagonal solver. Because this type of matrix inversion has the 
capability of simultaneously solving coupled systems of discretized (on a three point basis) equations, it has 
proven very useful in a number of (personally developed) problems. In addition to the implicit formulation 
described here, the block inversion method has been used as the base linear solver in a relatively complex, 
multispecies combustion modeling problem. This combustion problem exhibits strong coupling and non- 
linearity. Quasi-linearization and use of the block tridiagonal inversion method provided a workable solver, 
while decoupled scalar methods failed. For coupled marching problems, this seems to be an excellent matrix 
inversion method. 


3.6 Conclusions and recommendations for discretization methods 

This report has sought to describe the discretization using implicit finite difference methods and the 
numerical solution for a scalar, parabolic, linear partial differential equation. The areas described include: 

1. Derivation of a variable step, compact, implicit finite difference algorithm applied to a general, linear heat 
equation. Relationships to other methods are also considered. 

2. Resolution of singular behavior using analytical decomposition 

3. streamwise differencing. 

4. Block tridiagonal matrix solution algorithm. 

Many of the areas analyzed above were explored to try to develop methods that satisfied necessary constraints 
in the face of singular behavior. The linearity of this problem permitted the use of a relatively elegant and 
computationally successful decomposition to eliminate the singularity caused by the initial conditions. 

Since the ultimate intent of the underlying research has not been to make explore solution methods, but to 
provide a computational fluid dynamic tool, little numerical experimentation has been performed. Certainly, 
though, considerable experience was obtained in deriving and implementing these methods. Extending the 
table given by Fletcher (1984) using our experience, strengths and weaknesses of comparable methods are 
summarized below: 
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AttributeVMethod 

Explicit finite 
difference 

implicit finite 

difference 

(compact) 

Finite Element 

Spectral 

Trial Solution 

local 

“pseudo” global; 
Roach (1972) 

local 

global 

Ease of coding 

very good 

good 

good 

fair 

flexibility 
(geometric 
description etc.) 

good 

good 

very good 

fair 

accuracy per 
unknown 

fair 

good 

good 

very good 

computational 

efficiency 

good 

good 

good 

very good 

boundary 

conditions 

good, low order, 
fair, higher order 

very good 

good 

very good for 
simple conditions 

Matrix inversion 

very good, low 
order 

good, higher order 

very good, block 
inversion 

good 

very good, FFT 

Main strengths 

economy and ease 
of coding 

high order 
accuracy with 
minimal support 

flexible geometry 

high accuracy 

Main Weaknesses 

extension to higher 
order boundary 
conditions 

coupled 

formulation and 
inversion 

complexity, 
extension to non- 
linear problems 

Geometric 

inflexibility, 

complexity 


TABLE 3.2 Comparison of discretization methods. 


Table 3.2 provides a brief and perhaps somewhat prejudiced comparison of partial differential equation 
numerical solution methods. Newer methods, such as those based on wavelets, have not been included, since 
the author has very little knowledge of them. Other more classical methods, such as, traditional Galerkin 
methods, have not been directly included. Traditional Galerkin methods tend to be very poorly posed for 
simple polynomial (non-orthogona!) bases. Choosing a more suitable basis, such as, Chebyshev polynomials is 
leading us back to spectral and pseudo-spectral problems. 

With the limitations in the report in mind, the author concludes, that implicit or compact differencing 
methods offer an improvement in accuracy over simple explicit finite difference or finite element methods. 
They are far more useful than higher order explicit methods, i.e. five and seven point support methods, due to 
the superior treatment of boundary conditions. It is the authors opinion that implicit methods offer a method of 
intermediate complexity and utility between explicit finite element and finite difference methods and spectral 
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methods. Streamwise (temporal) differencing using the three point backward method provides a fully implicit 
method while preserving second order accuracy. The choice of lower order with respect to truncation error 
compared to stability restrictions seems to be appropriate. Numerical experiments described in Section 5 will 
also help justify the use of a lower order differencing scheme used in the streamwise direction, due to relatively 
gentle changes in the direction. 

Though significant recommendations for further work are made in Section 5, it is be useful to make some 
comments on area for further research concentrating on numerical methods here. Probably the biggest area of 
concern, is the development of compact formulas and methods which satisfy the conservation or control 
volume constraint. Equations (3.38) and (3.39) represent initial progress towards this work. Stability 
limitations for higher order methods, both streamwise and cross-stream, remain of interest. Exploration of 
methods based on Richardson extrapolation, predictor correctors or Gear’s methods; Ferziger (1981) may be of 
use for streamwise (temporal) differencing. Resolution of singularities, which was the driver behind much of 
the theoretical work presented here, remains critical. Finally, handling of nonlinearities in a manner that takes 
advantage of the higher order differencing, i.e. quasi-linearization with iteration, needs to be addressed. 


3. 7 Nonlinear equation solvers 


The majority of the closed form solutions developed previously are written implicitly in the form of a 
system of governing algebraic equations. As such, they often need to be inverted to yield appropriate variables. 
In this section the solvers used to invert these problems are briefly discussed. The methods employed: 


Problem 

Method 

Comments 

Use in DREA 

variable 

Scalar 

Bisection 

Linear convergence, robust 

rarely used 



Secant 

Superlinear Convergence, may 

Slipline static pressure 




fail to converge 

matching; pressure 

■1 




constrained entrainment 

■ 

System 

Gradient 

Linear convergence, robust 

rarely used 



Broyden 

Superlinear Convergence, may 

Fabri choke, critical back 




fail to converge. 

pressure 

Mi$. 





Pcri/Poi 


TABLE 3.3 Comparison of non-linear equation numerical solution methods. 


Although, Table 3.3 indicates, that several of the methods are “robust” while others may fail to converge, 
experience indicates that the faster converging methods (superlinear) are to be preferred. Convergence failure 
is more related to physics, i.e. no solution exists, than it is to numerical inversion. 

Since there has been little original development or modification to these solvers we refer the reader to the 
excellent references by Johnson and Riess (1982) for the scalar secant method and Burden and Faires (1993) 
for the Broyden method for descriptions of their development. 
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3.8 Variable grid Simpson's ride integration 


As the DREA suite of codes have been developed, the need to be able to obtain an integral average value 
of any quantity across the mixing duct has been apparent. A clear example of this need is the computation of 
the integral constraint (repeated here for convenience): 


!' 0 <f>(x,y)dy = i 


(3.38) 


The quadrature method chosen must: 

1. retain formal fourth order accuracy, to be consistent with the fourth order accurate difference 
formulation incorporated previously, and 

2. be valid on a variable grid. 

To satisfy these two requirements the most obvious choice is a variable step Newton-Cotes formula. One that 
will obtain the fourth order accuracy desired is a variable step Simpson rule quadrature formula. This type of 
formula is easily derived in principle, but is rather tedious to implement. As such, the development method is 
merely outlined here. We may start by defining the variable grid with the now known function f, as shown in 
Figure 3.3. 
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i 

f j + . 
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j+i 
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j 
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x j-i 



FIGURE 3.4 Variable grid definitions. 


To proceed, one next fits a Lagrange interpolating polynomial (Ferziger (1981)) through these three nodal 
points. This yields: 


m 


(*. - xX*,., - x) (x,_, - x) (x Jt ,-x) ^ , n , m 

A,(a ,. 1 + A,J" /w " V, ' MA-.+A)'" 


At this point it is possible to integrate, f(x) between Xj.i and Xj+i to obtain a rather cumbersome quadrature 
formula. 


\f{x)dx = /(a, , A y+1 ,fj ,f J+l ) (3-40) 
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It can be shown that for constant grid spacing, this equation reduces to the expected “1/3” Simpson rule 
(Ferziger (1981)). The net integral is then the summation of the intervals covering the range of integration. 

This formulation is poorly posed however, if the number of intervals chosen is odd. In this case it is 
formally necessary to go to the extent to derive an equivalent variable step 3/8 Simpson rule formula, which is 
even more tedious. The development is identical, except the quadrature stencil now involves four point 
support (and three panels), fit a Lagrangian cubic and proceed as before. A tedious procedure to implement 
but not conceptually difficult. 

Before concluding the formal discussion of integration methods, it is worth noting that since this 
methodology was based upon local polynomial curve fitting (i.e. Lagrange polynomials), the resulting 
quadrature formula is a very poor estimate in regions where the gradient is large. Indeed, this situation can and 
does occur in our flow field near the Heaviside step function inputs. As a bound on this error, consider the 
error associated with a 1/3 Simpson rule on a constant grid: 


/ (4) ( 7 /) 


Simpson 


90 


(3.41) 


Clearly, in regions where f 4) , i.e. the fourth derivative is large the error will be very large. 

In Section 5 it is shown that the integral constraint, as computed using our quadrature formula , is indeed a 
very poor estimate of the actual constraint near the singularity , for this reason. Clearly, the appropriate way to 
resolve this situation would be to introduce a combined numerical analytical quadrature formula. Fortunately, 
however, local integrations near the singularity are not of great interest. As such, it is acceptable to ignore this 
problem. This completes our discussion of the numerical solution methodologies. 
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4. EXPERIMENTAL METHODS 


In this section, inviscid, confined supersonic and subsonic stream interactions are studied using a 
hydraulic analog experiment. Specifically, the length of the first minimum location of the slipline is 
measured. Reasonable agreement is shown between a theoretical model and analog measurements. A 
method for extending the isentropic limited hydraulic analogy to model these inherently non-isentropic 
flows is also developed. This work is intended to aid in the general understanding of aerodynamic 
ejector/mixer nozzle operation, as well as, providing data for confined expansion compression surface 


4.1 Introduction 

Though free jets provide a useful limiting form of the ejector flow, the effect of confinement on the flow 
field is important. Thus measurements confined jet interactions provide data to validate the models developed 
in Section 2.3.2. Additionally, delineation of the inviscid flow field structure within an ejector nozzle has 
design implications for acoustic, structural, fluid mechanic reasons. 

In this portion of our study, a hydraulic analog/gas flow model of a 2-d ejector is used principally to locate 
and understand the structure of the inviscid multiple stream flow (principally the location of the aerodynamic 
choke). Figure 4.1 shows photographs of the hydraulic analog ejector model and slipline structure. As a basis 
for analytical comparison, formulas that describe wavelength of the slipline are developed for confined two- 
dimensional and axi-symmetric, supersonic-subsonic jets in De Chant et al. (1996) and Appendix C. 
Modifications have been made to this solution to extend its validity to include internal flows (De Chant et al. 
(1997)). This more general relationship provides internal flow field information for a supersonic ejector flow 
field, since the location of the first maximum expansion, approximately corresponds to the location of the 
aerodynamic or Fabri-choke, Fabri and Siestrunck (1958) or Addy et al. (1981). 



FIGURE 4.1 Water table analogy and expansion/compression structure. 
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In general, a formal mathematical analogy between two physical systems may be drawn whenever the 
non-dimensional governing equations and boundary conditions for the systems are equivalent. Examples for 
the Laplace equation include: electrostatic, membrane, Heleshaw and others (Bear (1972)). The inherent 
weaknesses of these methods are increased uncertainty due to additional theoretical basis. The analog 
discussed in this paper is a hydraulic analogy between gas flow and hydraulic flows. An analogy exists 
between thin, horizontal flow water under the influence of gravity and 2-d compressible, potential gas flows 
and this analogy has been utilized by a number of workers. Prieswerk (1940) performed some of the first 
practical studies using the analogy. Later workers include Shapiro (1954), Loh (1959) and (1969) and more 
recently, Pal and Bose (1993). 

Normal shock waves are approximately modeled by hydraulic jump phenomenon (Loh (1969)). The 
classical limitation for the hydraulic analogy, namely an effective specific heat ratio, y=2, is less restrictive if 
the analogy is used to verify an analytical tool for which it is possible to specify any convenient specific heat 
ratio. 


Gas Flow 

Hydraulic Flow 

Implication 

T >‘- r > 2 

ho . / 

— = 1 + — Fr 

To hp 

T 2 

d(pu) ( d(pv) 

h 2 

d(hu) | d(hv) 

T h ’ 7 

M 2 = Fr 2 ; ^ = ^ 
P h 

dx dy 

P 1 = P 1 To_ 

P P T 

dx dy J 

p UJ 


TABLE 4.1 Summary of hydraulic analogy. 


The key strength of analog methods is that they are less expensive relative to "full scale" measurements, in 
terms of, money, time and expertise. Table 4.1 is a summary of the relationships that may be drawn between 
the hydraulic and the analogous gas flow. 


4.2 Experimental facilities 

Practical application of the hydraulic analog is dependent upon the availability of a facilities where a thin 
hydraulic flows of water under the surface of gravity may be obtained, as well, as studied and measured. The 
Emil Buehler Aerodynamic Analog (EBAA) in the Mechanical Engineering Department at Texas A&M 
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University provides such a facility. Figure 4.2 provides a schematic of this facility. Detailed descriptions of 
the facility are provided by Felling (1992) and Lawton (1989). 


Water Table Wall 


— > 


Primary Reservoir 


Nozzle 


— > 


> 


— > 
— > 
— > 
— > 
— > 
— > 
— > 


Sharp Edge Weir 
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FIGURE 4.2 A schematic representation of the water table apparatus. 


The basic flow components include: 

• Pump 

• Primary and secondary reservoirs with flow rate measurement 

• Primary CD nozzle, and contoured secondary inlets 

• Ejector model test section 

• Return reservoir 

Originally, the EBAA facility was designed for single stream flows. Felling (1992). The EBAA facility 
was used with modifications to permit multiple stream flow rate measurement. A two stream manifold and gate 
valve array used in previous work was modified. The manifold, flow measurement, and piping array provided 
the basic source or reservoirs for the two streams. The large size of the EBAA reservoir provides very clean 
uniform flow for the secondary stream. The primary reservoir due to its smaller size has both a perforated “T” 
fitting and a series of three honeycomb screens to reduce primary reservoir turbulence. 

The EBAA facility was used with modifications to permit multiple stream flow rate measurement. A two 
stream manifold and gate valve array used in previous work was modified. These modifications included the 
installation of three Omega rotameter flow measurement devices. Specifications include: 

1 . Primary Stream: Single, 60 GPM Rotameter, Model FL75K, 2.5% full scale, 0.5% repeatability. 

2. Secondary Stream: Double 30 GPM Rotameters, Model FLD, 2.5% full scale, 0.5% 
repeatability. 


With primary and secondary reservoirs available, the two streams are mixed in a manner that simulates the 
ejector arrangement. Starting with the primary reservoir, a minimum length, converging-diverging nozzle was 
designed using the method of characteristics, see Shapiro (1953) and Anderson (1985), modified for hydraulic 
flows. The minimum length nozzle, Anderson (1985), accelerates the flow using a single “sharp ’ edged 
Prandtl-Meyer expansion. The exit flow is straightened using appropriately designed nozzle walls. These 
walls were designed using a graphical form of the method of characteristics method. Prandtl-Meyer expansion 
tables for heat capacity ratio y=2.0, required to satisfy the analogy, were generated and used in this 
computation. The nozzle throat was chosen to be 7.62 cm (3.0 in.) wide. For the design exit mach number, 
M=2.0, and throat, a nozzle contour was developed. The graphical value of the exit area 1 1.43 cm (4.5 in.) 
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compares well with the one dimensional value of 10.97 cm (4.25 in.) computed from the one-dimensional 
Mach number area relationship (isentropic mass conservation, y=2.0) 


1 throat 




2 1 , ' 
-0 + -(M e ) 2 ) 


(4.1) 


Comparing the graphical versus analytical exit area, the relative error here is approximately 6.0%. The 
subsonic (or subcritical) section of the flow was designed using large radius, smooth entrance inlets. The 
importance of large diameter inlets is emphasized by White (1986). The same type of inlet is used in the 
secondary flow field to minimize head losses. The two streams are brought together at the nozzle splitter 
plates, which are 0.3 1 75 cm (0. 125 in.) in thickness. The splitter plates induce local wakes, but these are small 
relative to the inviscid phenomenon of interest and may be neglected. 

The primary stream Mach/Froude number is a key parameter needed to model the inviscid interactions. It 
may be measured by employing the local Bernoulli relationship written between the primary reservoir and the 
primary nozzle exit plane: 


, , lu 2 

ho - h + ~ — 
2 8 


r. u n/ho ' 

Fr = 7 = [2(—-l)]i 

(gh)~2 h 


(4.2) 


This measurement requires measurement of the local water height As indicated in Shapiro (1954), these 
measurements may be made most conveniently using a needle probe arrangement. A simple, but effective, 
measurement tool consists of a beam (concrete level) with a pin probe (modified combination square, accuracy: 
±0.08 cm, ±0.0313 inches). More recently, Pat and Bose (1993) have used a flow visualization method to 
measure water depth. Though, this method has the advantage of being non-intrusive, it is not practical to apply 
to the two stream problem studied here due to significant three-dimensional disturbances inherent in a mixing 
flow. The pin probe introduces an additional degree of error, but remains the most practical and direct method 
for estimating local Froude numbers in the mixing flow. 

The location of the interface, may be readily observed by the marked change in water height at the 
interface. As seen from Figure 4.1, the interface is relatively well defined since it is caused by a rapid 
transition from supercritical flow (Fr> 1 ) to subcritical flow (Fr<l). Length was measured using the same beam 
and probe arrangement. Though the interface has small scale unsteadiness, it is of sufficiently small amplitude, 
less than 5% of the total inviscid structure length, and sufficiently high frequency, that an observer takes a 
“naturally” time averaged reading. 


4.3 Theoretical model and extensions of the hydraulic analogy 

This section summarizes the analytical model described in Section 2 and Appendix C, which is used to 
predict the geometry (length of first critical) of supersonic/subsonic confined and unconfined jets. These 
models provide a basis to understand the analog measurements. Then, the limitation of the hydraulic analogy 
in multiple stream flows is discussed and modifications to the analogy are developed using first order analytical 
arguments. 


4.3. 1 Interface slipline geometry (first critical location) 
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Following (De Chant et al. 1996) and Appendix C, relationships describing the slipline interface between a 
supersonic flow and subsonic flow free jet flow are merely quoted. The principle parameter is discussed is Xc/2 
which is the dimensionless location of the first minimum in the expansion/compression system resulting in the 
flow. See Figure 4.3. This value for both axi-symmetric and two dimensional flows is described by the 
equations below: 


free stream 
(2);U,M 2 <1 


Xc/2 first critical 
< 


> 


nozzle radius, y=l 


“ — . 

\\ " 

u 


>t 


^ interface slipline 


FIGURE 4.3 A schematic of the free jet problem geometry and definitions. 


•axi-symmetric: 

j = 1.22P, S -.22P„ - 1.22(1 + -j OP,.-22ft„ (4.3) 


•two-dimensional: 

f - 2p„ -P„ - 2(1 + jd) P„ - P„ < 44 > 

where p, =(M ls 2 -l) 1 ' 2 and ^ m =(Pu) 2 [l-U,/U ls ]. Equations (4.3) and (4.4) are the final modified solutions for 
the critical slipline location. These relationships with the lowest order terms, i.e. Pi s , U )s , based upon local 
static pressure matching for free jet problems and a modified method internal flow control volume theory, 
comprise our solution for the location of the first minimum of the slipline. By symmetry, half this distance is 
approximately the local maximum of the primary expansion and corresponds to the location of Fabn’s 
aerodynamic choke. Of course, only the 2-d formula is useful here since the analog is strictly valid for 2-d. 


4.3.2 Confined flow modifications 

Ejector/mixer nozzles are by their very nature confined or interna] flow problems as indicated by Figure 4.4. 
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FIGURE 4.4 Description of the intemal/ejector flow geometry and definitions. 


As such, modifications to the basic inviscid flow analyses of Section 2,3.1 are necessary. This modification 
takes the form of estimating the base flow or lowest order quantities in a different manner. Consider the 
expansion 


u d<j> 

— • = ; + —£•+ 

U is dx 




(4.5) 


where, U )s , is the velocity at the interface. In De Chant et al. (1996) the slipline velocity, U ls , was estimated 
using a static pressure matching methodology. To account for the internal flow effects in the current paper, it 
will be necessary to fall back on a more complete control volume methodology. Consider the two basic 
problems for an internal flow problem: 

(1) Large area ratio, unknowns: M ls , M 2s , Ui s 

(2) Small area ratio (Fabri problem), unknowns: Mi s , M 2 , U ts 

It is reasonable to expect, that two closures might be necessary to model large area ratio and small area ratio 
problems, in that, the limiting case of an infinite secondary flow field (free jet) cannot be analyzed at all using 
control volume methods. For a free jet, conditions are based solely upon local static pressure matching. This 
statement justifies the two modeling procedures: large area ratio problems are analyzed using static pressure 
and mass conservation, while for A 2 /Aj on the order of one, both mass and momentum conservation must be 
demanded. For a choked secondary stream, M 2s =l 5 this becomes Fabri’s problem. Analytical details of this 
work are provided by De Chant et al. (1997), and Section 2.3.2. 


4.3.3 Multiple stream extensions and limitations of hydraulic/gas flow analogy 

The hydraulic analogy as described in the introduction has been written for a single stream of water or gas. 
In this section, the analog is approximately extended to be capable of multiple stream flow simulation. This is 
done by first developing the classical analog with particular attention to restrictive assumptions. Then by 
examining a multiple stream flow the single stream limitation of the analog will be demonstrated. To 
overcome this limitation an approximate correction is introduced. The importance of this approximation is 
determined by comparison to experimental results. 

To understand the formal limitation of single stream, isentropic flow, the basic hydraulic/gas flow analogy 
is developed here. Consider the conservation equations for the two systems compared in Table 4.2. 


122 


Conservation Equation 

Hydraulic 

Aerodynamic 

mass 

V • ) = 0 

Po Uo 

v. (——) = 0 

ho u„ 

energy 


t ° = i + 1 t±m 2 

T 2 

momentum 


1 

ii 

-'fN- 

+ 


TABLE 4.2 Governing equations for hydraulic analog. 


Finally, although it has no counterpart in hydraulic flows, the state equation for the aerodynamic gas flow is 
written: 


p_ = _pL 

Po Po To 


(4.7) 


The relationships in Table 4.2 and equation (4.7) must be satisfied if an analogy is to exist between 
hydraulic and aerodynamic gas flows. By inspection, all these relationships may be simultaneously satisfied if 
one demands: 


p T h 

Po To ho 


(4.8) 


By either state or momentum it is necessary to write: 


£-rf/ 

Po h o 


(4.9) 


It should be noted, though, that this relationship is equivalent to invoking an isentropic assumption with y=2. It 
is apparent, that the isentropic assumption will limit the validity of the hydraulic analogy for any flow in which 
isentropic conditions are not well approximated. 

To emphasize this limitation, consider the case of two streams of different total conditions separated by a 
vortex sheet or slipline. By the definition, for a slipline to exist, the local static pressure is matched. Thus for 
two streams at a slipline, we may write: 


P,s = Po,[l + L ^-M, 2 JTr = p 2s = p 02 [l + M 2 2 


7_ 

i-r 


(4.10) 


Eliminating the static pressure and introducing the definition of total temperature yields: 
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(4.11) 


P 02 r Th -J- r P 02 fJL 

-f = [~^T Ir-if— Jh 

* oi 1 2 s l oi 


Now, alternatively, consider the total pressure definition required to satisfy the analogy: 


P02_ 

Pot 


j-T 02 -J- 

= hr 

1 oi 


(4.12) 


Hence, it is clear that except for the case of precisely matched static conditions (which is the trivial case of a 
single stream) that the hydraulic analog cannot precisely represent multiple stream flows. Our task is then to 
develop an extension to the hydraulic analogy that better represents these flows. 


4.3.4 An extended water table closure model 

The previous relationship, equation (4.12), indicates that the hydraulic analogy requires an isentropic 
relationship between the streams to be strictly valid. The definition of the total pressure at the interface 
described by equation (4.10) could be directly used, but it is in terms of the local interface static temperature 
ratio (and by analogy the local interface water height) which is difficult to measure. To understand this 
limitation from another point of view, a more physically based equation relating the total pressure jump to the 
total temperature using a local particular energy integral is developed. 

The closure relationship may be developed by considering Crocco’s law, White (1974), a particular energy 
integral, valid for small pressure gradient and Prandtl number equal to one: 


T01-T02 , 
dT 0 ~ du 

U1-U2 


and the definition of total temperature: 


dTo dT , du 

— - T +(r-i)M- 


Solving for the static temperature yields: 


dr - r , To V U, ' U2 ) n \ y ~ 1 \{>)1 dT ° 

' rp - L( rp rji )( )-(l + —r-M )] — 

1 1 01- 1 02 Uo 2 1 0 


The product term in the brackets may be simplified by assuming that: 


i+i* 
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(4.13) 


(4.14) 


(4.15) 


(4.16) 
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which is true for large variations in total temperature and interface velocity. These are precisely the type of 
flow conditions where the isentropic assumption needed by the analogy is at its weakest. Using a locally 
constant Mach number i.e. average between streams permits integration of equation (4. 1 5) to yield: 


III = (hi 

T 2 s To/ 

Note that for incompressible flow M ave =0, that the proper limit is obtained: 

JjlJ hi]' = hi 

T 2s ^7V Toi 


(4.17) 


(4.18) 


Hence, the approximate relationship for the total pressure in terms of the stream total temperatures valid for 
multiple stream flows is obtained: 


hi 

Pot 


To J 


Simplifying the exponent of equation (4. 1 9), yields: 


y f(r-i) 
y — 1 . 2 


M 



(4.19) 


(4.20) 


The classical analogy, though, requires that equation (4.20) take the constant value of two. From equation 
(4.19) it is clear that satisfaction of analogy for multiple stream flows is Mach (or Froude) number dependent. 
It is our conclusion that in air, where y= 1.4, that only for an interface Mach number of: 


M„ 



*1.69 


(4.21) 


can the analogy be expected to hold for non-isentropic flows. Additionally, since it has been shown that for 
multiple stream flows, the exponent in the pressure relationship is Mach number dependent, slight modification 
of the exponent in the classical closure i.e. 2— >2±e, with e a small correction term, may be expected to extend 
the validity of the model. To satisfy the other relationships in the analogy, y=2.0, should be used for all 
computations. This would yield an average Mach number value: 

M *(-r * 1.4142 (4-22) 

* ave l 2/ 
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Ultimately, as long as, either equation (4.19) or (4.20) is approximately satisfied, the analogy can be expected 
to provide reasonable results even for multiple stream interactions. 

Later in Section 4.4, the validity of this correction strategy is ascertained by comparison of the analog data 
for a case where Mi=2 and the secondary stream chokes. For this flow, the correction factor, e, may be 
estimated. Neglecting the acceleration of the primary stream, but recognizing that the secondary stream is 
sonic yields a local average Mach/Froude number, M avc ^(2.0+1.0)/2=L5. To approximately satisfy equation 
(4.22) and, more importantly match the estimated local Mach/Froude number, the classical exponent value of 
two is modified 2->2±e, with e=l/8. From equation (4.22) this would yield an average Mach number 
M avc =[(2.125) 2 /2] a5 = 1.5026. As will be seen, this correction term provides a slight improvement to the model 
for multiple stream flows. A summary of the basis and form of this correction is provided in Table 4.3. 


Gas Flow 

Hydraulic Analog 

Extended Analog (M ave =l .5) 

Pot \T 2 ) VTflyJ 

Po2_(ho2^\ 
Po,\h 0 , J Y ~ 

*-r*ar tm i 

Pot v/j0K 8 


TABLE 4.3 Extended hydraulic analogy. 


4.3.5 Error and error propagation for the hydraulic analogy 

The previously developed analyses have essentially assumed that the flow of water in the hydraulic 
analogy is a 1-d, frictionless flow. Three dimensionality and frictional losses are significant for the two stream 
mixing problem described here. As a way to understand the effect that losses have upon the analogy consider 
the Froude (or by analogy Mach) number computation: 


Fr 


2 (^- 1 ) 


1 + f 


friction 


J 


(4.23) 


As shown in equation (4.23), friction has the effect of reducing the actual Mach number. Off setting this 
reduction is the presence of local three-dimensional effects and turbulence caused by the rapid acceleration of 
the flow through the nozzle. Biasing is also possible in the nozzle design or fabrication. Measurements of the 
exit Mach number for 20 separate cases using the pin probe and beam arrangement and equation (3) yield an 
average exit Mach/Froude number of 2,09, a relative error of 4.5%. 

An additional source of uncertainty comes from the actual measurement of the Froude number. As a way 
to quantify this error and its propagation, consider the uncertainty analysis method of Kline and McClintock 
(1953). For our problem, it is sufficient to consider the Froude or Mach number squared relationship: 

Fr ' m ifi~ i 2(! k' i)] (4 - 24) 
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Now, since both h 0 and h are directly measured, using the pin probe apparatus, the partial derivatives of the two 
measured quantities, h 0 and h are computed: 


cFr 1 2 dFr 2 

h ch \h)h 


(4.25) 


Introducing the ratios Ah</ho and Aho/h 0 where the Ah is the error associated with the pin probe measurement, 
squaring, summing and taking the root yields the propagated uncertainty: 


hj h 0 h 


Ah 


(4.26) 


The errors associated with the two measurements for a representative case are: 

• primary reservoir value: Ah</ho = 0.03 l(inches)/3.0 inches=0.01 

• secondary reservoir value Ah/h=0.03 1 inches/1 .0 inches=0.03 1 

Substituting these values into equation (4.26) yields a value for the Kline-McClintock uncertainty of 
approximately o) Fr 2 =20%. This value gives a rather pessimistic estimate of the error associated with the 
measurement of the Mach number in the primary stream. The actual error is probably considerably smaller, as 
indicated by the small relative error 4.5% shown for the primary Mach number. 

The good comparison between aerodynamic theory and the hydraulic analogy experiment demonstrate that 
the hydraulic analogy may provide a viable method for making quantitative measurements in support of 
numerical and analytical modeling efforts. Additionally, the rather generous error bound from the Kline- 
McClintock theory is shown to be just that, a pessimistic upper bound on the error. 


4.4 Results ; analog versus theoretical results and repeatability 

The fundamental result for this problem is the measurement of the critical location, Xc, of the slipline. 
Results for two cases are presented in Figures 4.5 and 4.6. The relevant operating conditions for these cases 
are summarized in Table 4.4: 


case 

Primary Mach No., Mi 

Secondary Mach No. , M 2 

Area Ratio: A 2 /Aj 

1 

2.0 

negligible 

9.60 

2 

2.0 

negligible 

3.235 


TABLE 4.4 Operating conditions for slipline geometry measurement. 
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Total Temperature Ratio T02/T01 

FIGURE 4.5 Comparison of experimental and model critical interface lengths, for A 2 /A,=9.60. 

Since both of the area ratios for the problems measured are relatively large, the closure model chosen is the 
large area ratio static pressure matching method, i.e. equation (2.40). Analog measurements have been 
consistently reduced using the extended model. Both the unmodified model and the extended model have been 
used to provide total pressure ratio input values into the control volume method. Additionally, a series of 
repeatability test were run, using dimensional water height combinations to yield the same reservoir ratio’s , to 
assess the effect of actual water heights upon the measurements. A significant, but acceptable, amount of error 
is introduced by varying this ratio, especially for significant differences in reservoir height, i.e. T 02 /T 01 small. 
For these cases, three-dimensional effects are very pronounced. Though not presented here, repeatability using 
similar reservoir ratios, but rerunning the overall experiment, is excellent. 
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FIGURE 4.6 Comparison of experimental and model critical interface lengths, for A 2 /A]-3.235. 


A series of general conclusions may be drawn from the Figures 4.5 and 4.6. 

• The theoretical model, De Chant et al. (1996) and Section 2.3.1., provides reasonable agreement with the 
analog measurements. This good agreement gives confidence in the simple model described previously. 
Alternatively, the analog provides a reasonable simulation of inviscid gas flow phenomenon including the 
choking phenomenon characteristic of supersonic nozzle flow. 

• The hydraulic analog provides an inexpensive, adequately repeatable model for studying inviscid, multiple 
stream flows. 

The simple correction model, i.e. modification of the total pressure ratio term 2— >2±s, with e=l/8, provides 
moderate improvement for the analytical model. 

4.5 Experimental conclusions 

Inviscid, confined supersonic and subsonic stream interactions have been studied using an 
experimental hydraulic analog model, as well as, an approximate, perturbation theoretical model. 
Specifically, the first critical location of the slipline is considered. Reasonable agreement is shown between 
the theoretical model and analog measurements, hence helping to validate the concepts developed in 
Section 2. Additionally, a method for extending the isentropic analogy to model these inherently non- 
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isentropic flows was discussed. This work was of considerable use in extending our understanding of 
aerodynamic mixer ejector nozzle operation. 
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5. RESULTS AND CONCLUSIONS 


5. 1 Results 

This section compares the computational methods developed previously to experimental 
measurements. Given the wide range of approximations inherent in this modeling effort, this type of 
comparison is absolutely mandatory. Since many of the supporting analyses have been partially validated 
along the way, the overall computational code, DREA is compared in detail. Included in these results are: 

• flow conservation and grid independence studies 

• integral quantities: pumping/entrainment and CFG 

• local quantities: streamwise and cross-stream profiles of Mach number, velocity, and static 
pressure. Additionally, we present limited contour information to get a view of the overall ejector 
phenomenon. 

It should be pointed out, that results are virtually always computed and presented in dimensionless 
form. The input files use dimensional quantities only to ease interfacing with existing codes. Spatial position, 
x and y are scaled by the inlet centerline to shroud distance: 


X non-dtim i / rw ^ non -dim 

m 


y 

m 


(5.1) 


where h(0) is the inlet shroud to centerline distance as shown in Figure 5.1. 


shroud 



FIGURE 5.1 Definition of non-dimensionalization quantities. 


Additionally, the dependent variables, both conservative and primitive variables, are scaled by the inlet mixing 
plane, area averaged conditions. These conditions are chosen, since they are readily available. In other words: 


fnon-4m( X >y ) = 


f(x,y) f(x,y ) 


fme + flO^. 


'(■^10 ^20 ) 


(5.2) 


20 


Note that the “non-dim” subscript is virtually always suppressed. 


5. 1. 1 Flow conservation and grid independence studies 
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In this section, some computational aspects of our simulation are considered. Here, we are principally 
concerned with the numerical accuracy of our simulation. To ascertain numerical errors, concentration is 
placed upon, the momentum flux, G=pu 2 +p (which is conservative) for a simple subsonic flow problem with a 
constant area shroud. Our interest is in the degree to which the integral constraint will be satisfied: 

\' 0 <p(x,y)dy = l (5.3) 

where, (p=G. The degree to which the constraint is satisfied is dependent upon several factors that are presented 
in Table 5.1: 


solution portion 

dependency 

comments 

numerical solution 

Ay (jmax) and Ax (i max ) 

truncation errors in streamwise and 
cross-stream direction 

analytical solution 

^max 

number of terms in Green’s 
function expansion 

numerical integration 

Ay (jmax) 

cannot resolve near field singularity 

shroud geometry 

number of wall points and length 

second order 


TABLE 5.1 Sources of numerical and truncation error for integral constraint. 


We start by considering a subsonic, constant area shroud problem the size of the constraint, again which must 
equal one, for a basic jmax = 20 and i max = 30. As is clearly shown in Figure 5.2 (and accentuated by the 
logarithmic scales), the near field value of the constraint, is poorly predicted. 
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FIGURE 5.2 Values of the integral constraint at various streamwise locations. 


This failure is directly related to the inability of the polynomial based, variable step integration 
routine, which is incapable of resolving the near field singularity. This limitation is exactly analogous to the 
inability to resolve the near field “Heaviside” step function using solely numerical means. This limitation is 
somewhat irrelevant, in that in a very short distance, the solution recovers the desired constraint value. Further, 
it should be emphasized, that the quadrature method is in error; not the underlying functional valu es at the 
knots, since these values were computed using our combined analytical/numerical method which is capable of 
resolving the singular behavior. 

Due to the fact, that the integration routine cannot be used in the near field it would be more sensible 
to consider the degree of error, as a function of grid parameters, in the constraint at the exit of the flow. 
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FIGURE 5.3 Relative error in satisfaction of integral constraint versus number of cross stream grid points 
evaluated at ejector exit plane for a subsonic ejector. 


As can be seen from Figure 5.3 the expected grid dependence is observed. The numerical error, does indeed 
does drive to a maximum (where round off issues may become important, more likely, though, the definition of 
the shroud geometry which uses a cubic spline basis is the key limitation). Additionally, a simple regression of 
the linear (on a logarithmic plot) portion of this plot. Figure 5.4: 


1E-2i 

'<■ • 

1E-3-- 

Relative error IE-4 - ' * 

IE-5 

1 E-6 J 

0 1E1 2E1 3E 

Number of cross-stream grid points 

FIGURE 5.4 A regression fit of relative error versus number of cross-stream grid points at the exit plane 
which demonstrates 4th order accuracy for a subsonic ejector. 
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clearly indicates the fourth order (exponent=-4.77; actually almost fifth order) dependency of the constraint 
upon the number of cross stream grid points. It is shown, that the numerical portion of the code is contributing 
very little error. Most of the error in this problems is being generated due to the integration routine itself. 



dimensionless streamwise location, x 
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0.8 


FIGURE 5.5 Contour plot of momentum conservation term, G, demonstrating rapid cross stream variation 
and gentle streamwise variation for a subsonic constant area ejector. 


Performing the same type of computations, except now holding the cross-stream coordinate constant 
and considering the effect of the streamwise coordinate shows the insensitivity of the modeling to streamwise 
gradients. This is of little surprise, in that, the most significant changes occur in the cross stream direction. 
Recall that it was indicated in Section 3 that this was the basis for accepting a method which is 4th order in the 
cross-stream and only 2nd order in the streamwise direction. Additionally, the numerical portion of the solution 
is virtually negligible for this problem. Variation in the streamwise direction is scaled at 0(1). As expected, 
and indeed found, the solver is virtually independent upon the number of streamwise, I mlx steps. The gentle 
streamwise change is shown in the Figure 5.5: 

Similarly, it is found that the total solution is virtually independent of the number of terms of the series 
used to evaluate the numerical portion of the total solution for G, i.e. N,*,. As shown in the pre vious figure, the 
relatively small size of the primary stream denotes a small value for h s , hence reducing the required value for n. 
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In fact, it is expected the numerical contribution for this case to be negligible. 

We do not desire to leave the reader with impression, that the numerical solver is unimportant. If one 
considers an alternative supersonic constant area shroud test problem and print out the values for (G m - 
GnunO/Gm, it is possible to see that, indeed the numerical contribution to the solution is non-negligible. See 
Figure 5.6. Additionally, the local streamwise changes in the number of terms used in the numerical solution 
are apparent, i.e. x=3. Further, as expected, the numerical contribution near the shear layer is small and 
increases upon approach to the boundaries (especially the y=l boundary). 
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FIGURE 5.6 Contour plot of (G^G^yG^ for a constant area supersonic ejector indicating non-trivial 

numerical contribution far from centerline; y=0. 


Performing the same types of grid independence tests for this problem, Figure 5.7: 
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FIGURE 5.7 Relative error in satisfaction of the integral constraint versus number of cross stream grid 
points evaluated at the ejector exit plane for a supersonic ejector. 


The constraint error for the supersonic problem is also decaying at a rapid rate. Regression and Figure 5.8 
indicate that the exponent for this problem is approximately, -5.5. 
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FIGURE 5.8 A regression fit of relative error versus the number of cross-stream grid points at the exit 
plane which demonstrates 4th order accuracy for a supersonic ejector. 
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Value indicates, that the numerical solver and the numerical integration algorithm are indeed yield fourth order 
accuracy (or better!). 


This completes our discussion of grid independence and conservative integral constraints. It is 
reasonable to conclude that the combined numerical/analytical solution is providing a consistently high 
accuracy solution. Of course, this conclusion must be tempered, by the reality that though the equations are 
being solved accurately and consistently, the value of the solution is only as good as it compares to physical 
solutions. This demands a comparison with experimental results (which may contain measurement and 
reduction error) but have the greatest chance of mimicking the physical prototype, in that they are 
measurements of physical processes and quantities. Comparison with experiment is discussed in the next 
section. 


5. 1.2 Summary of comparisons 

Here a family of test cases of increasing complexity is presented. The latter two test cases are from 
current industry/NASA sources. Comparison to this data is critical since it represents the forefront of 
aerodynamic ejector technology, which is the family of problems that DREA has been designed to model. 
Hence use of this data indicates that DREA can make a significant contribution to actual modeling problems. 

The list of ejectors modeled is presented in Table 5.2. 
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5. 1.3 Gilbert and Hill (1973) 


This is a two-dimensional, (planar splitter plate), subsonic (choked primary), variable area shroud 
ejector. Results are provided for pumping in Table 5.3 and streamwise mixing (including pressure, velocity and 
temperature fields). 


Entrainment 

experimental measurement 

DREA computation 

relative error 

w 2 /w. 

3.76 

3.69 

1.86% 


TABLE 5.3 Comparison of DREA simulation entrainment prediction with Gilbert and Hill (1973). 


The comment may be made that the entrainment prediction is very good. This type of accuracy is probably 
fortunate, rather, than something to be relied upon. Additionally, since the ejector is completely subsonic, it 
should be expected that the entrainment will be insensitive to the length of the shroud. Of course, this is an 
approximation. 

To gain a feeling for the flow field mixing, consider the normalized velocity, static pressure and static 
temperature. They are presented in Figures 5.9, 5.10 and 5.11, respectively: 


140 














FIGURE 5.9 A comparison between the DREA simulation and experimental data of Gilbert and Hill (1973) 
showing centerline velocity versus streamwise location with uncertainties estimated from the literature. 


We note that the simulation provides a reasonable comparison to the experimental data over the entire 
range. Some deviation is observed in the mid region of the flow, where variable area shroud effects are 
especially important (and not precisely modeled). Fortunately, the Gilbert and Hill (1973) report, relative to 
other data sets, provides information on the experimental error and, more importantly, information on the 
associated uncertainty. Gilbert and Hill estimate their uncertainty to be on the order of 4 %. This uncertainty 
estimate is the experimental workers best estimate of error propagation, as well as, an estimate of the degree of 
scatter from repeatability tests. Four- percent error bars are introduced into Figure 5.9. 
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FIGURE 5.10 A comparison between the DREA simulation and experimental data of Gilbert and Hill 
(1973) showing the centerline static temperature versus streamwise location. 


Comparison with the static temperature, shown in Figure 5.10 is adequate. Though the plot tends to “magnify” 
the error, the relative error is no more than 8%. 
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FIGURE 5.1 1 A comparison between the DREA simulation and experimental data of Gilbert and Hill 
(1973) showing wall static pressure versus streamwise location. 

Wall pressure values are well modeled as shown in Figure 5.11. Here we have introduced the low Mach 
number wall correction described in Section 2. Next a supersonic/subsonic mixing layer problem is considered. 


5. 1.4 Goebel and Dutton (1991) 

This is a supersonic subsonic mixing problem. The primary Mach number, Mi=2.5 and the secondary 
is M 2 =0.3. The secondary static conditions are controlled so as to force matching (or closely approximate it) of 
the static pressures at the interface. As such, no independent entrainment information is available. 
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FIGURE 5.12 Comparison between DREA simulation and experimental data of Goebel and Dutton (1991); 
dimensionless velocity difference profile, (u-u 2 )/(u r u 2 ); x=50 mm. 


It is possible, though, to assess the degree of mixing in the flow. Here, the interest is the width of the mixing 
zone. Comparisons are provided for two streamwise locations x=50 mm and 100 mm, Figures 5.12 and 5.13 
respectively, measured from the splitter plate. 
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FIGURE 5.13 Comparison between DREA simulation and experimental data of Goebel and Dutton (1991); 
dimensionless velocity difference profile, (u-u 2 )/(u r u 2) ; x=100 mm. 


As can be seen, the mixing width (slope of the mixing layer) is predicted very well. This is 
encouraging, indicating that, at least locally, the turbulence model is providing a correct estimate of the mixing 
rate. Unfortunately the location of the mixing layer is not as good. Several explanations for this minor 
disagreement are possible. One possibility is that there are wall boundary layer blockage effects, distorting the 
streamlines. More likely, there are small pressure imbalances (even locally) which result in expansion or 
contraction of the dividing slipline. 


5. L5 Fernado and Menon (1993) 

The next problem, considered is a supersonic (Mach 2.5 and, choked flow (Mach 1.0) mixing problem 
studied by Fernando and Menon (1993). This study involves a 2-d mixing layer shown schematically in Figure 
5.14. 
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locally constant area mixing 


Primary Mach=2.5 
secondary Mach=1.0 



FIGURE 5.14 A diagram of the basic experimental apparatus for Femado and Menon (1993). 


Due to the large change in shroud geometry, this flow is modeled as a straight wall (and our comparison is 
confined to early in the flow field). Figure 5. 15 shows a 2-d slot ejector, which is essentially a reference case. 
Comparison with experiment is adequate. Again, the physical location of the shear layer is not exactly 
predicted. Here, wave expansion/compression effects around the step are unavoidable, potentially modifying 
the interface slipline. Results are presented in Figure 5.15. 



FIGURE 5.15 A comparison between the DREA simulation and experimental data of Fernando and Menon 
(1993) showing the Mach number profile for a slot mixer. 
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The second case considered, is a vortically mixed convoluted, splitter plate. The three vortical parameters for 
this case are presented in Table 5.4. 


Vortical Reynolds Number 

Lobe total height/lobe wavelength 

Lobe total height 

Re V oit = 0*087 

h<A=0.428 

ho/h(0)=0.0061 


TABLE 5.4 Vortical parameters for Fernando and Menon (1993). 


Since the lobe penetration (height) is relatively small, the effect on overall mixing is slight. Locally, however 
the flow is locally accelerated near the chute causing the experimental jump in Mach number. 


2.80 


2.40 


JJ 


o 


2.00 


1 .60 — 


Fernando and Menon Experiment (1993) 
Conditions: Ml=2.5, M2=1.0; 
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FIGURE 5.16 A comparison between the DREA simulation and experimental data of Fernando and Menon 
(1993) showing the Mach number profile for a vortical mixer. 


From figure 5.16, it is clear that the computational model, though, cannot predict this local acceleration as the 
flow comes off the chutes (recall, that our model is averaged in the z-direction). Alternatively stated, the 
computational model cannot detect local primary or secondary chute disturbances. The model does, however, 


147 





include the enhanced mixing effect upon the average stream conditions. Considering the overall Mach number, 
though, we see that the mixing zone width is well predicted. 


5. L 6 Arney and Lidstone 

The next family of predicted problems involves advanced nozzles tested at Boeing. These are combined 
supersonic primary and subsonic secondary stream mixing flows, Mi/M 2=2.86. A large family of tests has 
been considered here. Concentration is placed upon several of them. The Boeing nozzle tests are a simple, 
slot ejector problem. The intention of these tests though is provide local information, since in close proximity 
to a lobe , the flow acts like a 2-d shear layer. 

Though it is not totally clear from the report, it seems apparent that secondary quantities have been 
controlled, hence secondary entrainment is not an independent code predicted quantity. These cases, though, 
do provide an excellent test of mixing within the duct. Experimental data was taken using LDV Laser-Doppler 
velocimetry. The conditions simulated in this report are listed in Table 5.5. 


Figure 

Total Temperature ratio; 
T 02/To i 

Total Pressure ratio; P02/P01 

Streamwise Location, 
x/h(0) 

5.17 

1.0 

3.5 

1.05 

5.18 

1.0 

3.5 


5.19 

0.49 

3.5 

1.05 

5.20 

0.49 

3.5 

4.2 

5.21 

0.332 

3.18 



TABLE 5.5 Parameter summary Amey and Lidstone. 


Additionally, limited experimental error estimates have been provided, especially for the LDV system. No data 
uncertainty estimates have been provided at all. 

The first set of comparisons Figures 5.17 and 5.18 are for the conditions with NPR (nozzle pressure ratio, 
i.e. poi/p«) of 3.5. Total temperatures are matched for this case. 
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FIGURE 5.17 A comparison between the DREA simulation and experimental data of Amey and Lidstone 
using a dimensionless velocity profile for a slot mixer, with x/h(0)=l .05. 
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FIGURE 5.18 A comparison between the DREA simulation and experimental data of Amey and Lidstone 
using a dimensionless velocity profile, a slot mixer, x/h(0)=4.2 and matched total temperatures. 


Alternatively, a hot flow set of measurements (implying dissimilar total temperatures), Figures 5.19 and 5.20, is 
considered: 
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FIGURE 5.19 Comparison between DREA simulation and experimental data of Amey and Lidstone; 
dimensionless velocity profile, slot mixer, x/h(0)=l .05. Hot flow case; To2/Toi = 0.4890. 


and at a location farther downstream: 
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FIGURE 5.20 A comparison between the DREA simulation and experimental data of Amey and Lidstone 
using a dimensionless velocity profile, slot mixer, x/h(0)=4.2 and To2/Toi=0.4890. 


Finally a further set of hot flow measurements, i.e. strongly dissimilar total temperatures, are shown in Figure 
5.21: 
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0.00 0.40 0.80 1.20 

Cross-Stream Location; y/H 

FIGURE 5.21 A comparison between the DREA simulation and experimental data of Amey and Lidstone 
using a dimensionless velocity profile, slot mixer, x/h(0)=4.2 and T 0 2 /Toi=0.33 17. 

Equally important are thrust considerations. In the following figure. Figure 5.22, CFG (Gross thrust 
coefficient) and pumping (entrainment) are provided. 
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FIGURE 5.22 A comparison between the DREA simulation and experimental data of Amey and Lidstone 
showing entrainment and gross thrust coefficient for a slot mixer. This is a hot flow case. 
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While it is clear that the trend is correctly predicted, the values are in moderate to large relative error. Of 
special importance (this is related to the flight regime of interest) is the mid location of these graphs (i.e. NPR 
2-4), where one may note relative errors on the order of 11-23%. The normalization values used are the CFG 
and entrainment values evaluated at an NPR of 3.0, respectively. 

The source of this error is somewhat hard to estimate. Since the velocity field has been matched 
adequately as presented in Figures 5.17-5.21, only the density computation or (and we believe this to be most 
likely error source) the inlet cross-sectional areas are incorrect. Recall that, since the secondary stream was 
controlled by upstream valves, the control volume pumping routine has not been employed, and is therefore not 
the cause of this discrepancy. That some form of inconsistency may be contained within the data is also 
supported by the CFD simulations performed by R. K Henke (1994) who showed good velocity match, but an 
unrealistic expansion\compression interface structure (the static between stream is matched, hence the “shock 
diamond” system is erroneous). Henke notes the problem, but offers no explanation. Our computations 
indicate that the static pressure is indeed matched. The possibility that boundary layer blockage is biasing the 
computations must also be considered. 


5. 1 . 7 Thayer et al. 

The final set of cases are vortical supersonic/subsonic (subsonic at exit) cases. Figure 23 compares the 
CFG and entrainment computed by DREA with that from the test data. 
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FIGURE 5.23 A comparison between the DREA simulation and experimental data of Thayer et al. with 
estimated experimental error for a vortical mixer. 

As can be seen from Figure 5.23, the error here is relatively small for the CFG, 10-12%, but ranges 
from negligible to 100% for entrainment. The temperature corrections consist of multiplying CFG and 


156 





entrainment by (T 02 /T 0 i) 1/2 , respectively. The normalization values are chosen at an NPR of 3.5. Of special 
interest is the middle range which again shows tolerable error. The error bars represent a 0.25% 
measurement error. No uncertainty estimates are available. 
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FIGURE 5.24 Normalized Mach number field for vortical mixer, chute acceleration neglected. 


It should be noted that flow field geometry fin particular ejector width) was not supplied and had to be 
estimated, thus introducing potentially large entrainment and CFG error. The normalized Mach number 
field, M(x,y)/M|, is presented for this case in Figure 5.24. 

These tests are for a vortical mixer ejector nozzle with the mixer parameters listed in Table 5.6. 


Vortical Reynolds Number 

Lobe total height 

Re vort =2.0 

h</h(0)=1.0 


TABLE 5.6 The vortical parameters for the Thayer et al. experiment. 


This is clearly a stronger mixer than described previously. Recall that the DREA simulation cannot 
predict local Mach number acceleration of the chutes. As such internal static quantities are poorly predicted. 
To understand this limitation more completely, consider the dimensionless mixing duct static pressure for this 
problem, with the modification of 16% higher inlet Mach number, shown in Figure 5.24. 
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0.00 1.00 2.00 3.00 
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FIGURE 5.25 A comparison between the DREA simulation and experimental data of Thayer et al. using a 
dimensionless static pressure p^/p* for a vortical mixer. 

Figure 5.25 indicates that the local static pressure near the inlet is not correctly predicted. A higher inlet 
Mach number would permit better pressure matching, but at the expense of even worse prediction of the 
integral quantities. Notice also, the weak shock near the exit of the ejector and the diffusion of the high 
pressure flow from this point. This is illustrated by the dimensionless contour plot (here, p(x,y)/poo) in Figure 
5.26: 
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FIGURE 5.26 Contour plot of p(x,y)/p«o, vortical mixer. 


Notice the high primary stream terminating in a long shock train. As expected, the static pressure reaches a 
maximum at the sonic point (location of the weak shock system. This correspondence is shown in Figure 5.27 
by the normalized Mach number, M(x,y)/M l0 , for this case: 
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FIGURE 5.27 Contour plot of normalized Mach number, M(x,y)/Mi 0 , vortical mixer. 


Recalling that the inlet flow characteristics have been modified, i.e. local increase in average Mach number, it 
should be expected that a poor match to the integral quantities, entrainment and CFG, would exist. This is 
indeed the case as shown in Figure 5.28: 
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FIGURE 5.28 A comparison between the DREA simulation and experimental data of Thayer et al. showing 
entrainment and CFG, using a vortical mixer and modified initial conditions to account for mixer lobe 

acceleration. 
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Once again, it is apparent that some error has been introduced. Clearly, the entrainment is suppressed for the 
higher Mach number case. The failure of the code to predict the local acceleration due to the lobes of the 
mixer has caused most of this discrepancy. Methods to extend our model to begin to account for these effects 
are discussed in Section 5.2. 


5.1.8 NASA Lewis design study , 1996 

As a support for design trade studies and as a demonstration of another solution of the code, preliminary 
results from a NASA “in house” parametric study have also been included. This study has been chosen, in 
spite of the fact that no experimental data is available, since it illustrates an ejector running in an aerodynamic 
or Fabri choked mode. Presented in Figure 5.29 is the normalized flow field Mach number: ' 

Figure 5.29 shows a supersonic stream initially mixing with a subsonic stream. At x=1.6 a rapid 
acceleration of the secondary stream is occurs. This is the estimated (using the interface streamline model) 
location of the aerodynamic or Fabri choke. It should be noted that the rapid acceleration is an artifact of our 
analysis and should indeed be more gentle. This is a clear indication of the limitations caused by being unable 
to model the wave phenomenon (expansion/compression system) within this flow. The blending function 
approach discussed in Section 2 was found to produce unrealistic results. 



FIGURE 5.29 Contour plot of normalized Mach number, M(x,y)/Mi, for Fabri choked flow, NASA Lewis 

design study. 
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5.1.9. Preliminary computational time studies 


The NASA design study is also useful, in that it is potentially the most complex type of flow problem we 
are likely to run. As such some crude computational time studies have been performed on a Pentium 90 PC 
platform using Microsoft FORTRAN V. 1.0 (all compiler optimization options are set for default). These are 
user time studies for code, i.e. time experienced by the user from initial submission of the job to completion 
and are probably I.O. limited rather than CPU limited. Results are presented in Table 5.7: 


Streamwise grid 1^ 

Cross-stream grid J max 

Time (seconds) 

0 

0 

3.15 

20 

10 

4.66 

20 

20 

5.12 

20 

30 

6.18 

30 

10 

5.13 

30 

20 

6.5 

30 

30 

10.1 

100 

30 

20.0 


TABLE 5.7 Elementary computational time studies. 


We note that the first case represents a run computing integral quantities, such as entrainment and CFG only. 
To our satisfaction, we note that the results of this computation compare favorably to the time required by Fung 
(1995) of about 4 seconds on an IBM RS6000 workstation, to obtain their “momentum mixedness” parameter 
and integral quantities. At the cost of on the order 0.5 to 1 .5 seconds more, our analysis provides flow field 
information at 200-400 grid points within the shroud using a less powerful (and vastly more accessible) 
computational platform. Indeed, though comparison of computational platforms is difficult (for example 
integer performance versus floating point performance) it is reasonable to expect that a IBM RS6000 560 
platform should be on the order of 3-4 times as fast for floating point computations. With this in mind, it is 
apparent that the DREA code is providing vastly more information a t essentially half the cost ] 

For comparison purposes it would be useful to have an empirical formula relating grid points to 
computational time. Since I max is not a strong variable in this determination it is held constant, i.e. I max =20 and 
we compare J max to computational time as shown in Figure 5.30: 
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FIGURE 5.30 A regression fit of the computational time on PC Pentium 90 platform versus the number of 

cross-stream grid points. 


Clearly, the cost of larger numbers of computational points is found to reasonable and increases in a linear 
manner. Of course, this empirical relationship could have well been fitted by a power law expression: 
y=3.3 lx 1 ' 02 , but again it is apparent that this is a virtually linear relationship. Additionally we can see that for no 
more twice the cost of a control volume solver, which provides no length scale or internal mixing information, 
we can use this model to obtain approximately 600 points of information within the mixing duct! 

At this point our discussion of the validation cases is concluded. These six cases and the parametric case 
have helped us ascertain that the code seems to provide reasonably accurate mixing or flowfield, entrainment 
and thrust information. We have also observed the limitations inherent to the model. Both of these topics are 
developed in the next section. 


5,2 Conclusions and recommendations 

The mathematics, numerics, supporting experimental research and comparisons of the DREA simulation 
code have been presented. It is now important to draw conclusions for this research, as well as rigorously 
define extensions and areas of further work to improve or supersede this work. 

5.2,1 Conclusions 

In this section the success/or lack of success achieved by this development project is discussed. There are 
two major components to this research. 

1. General goal : To develop a new and useful family of approximate fluid dynamic equations which provide 
alternatives to classical hierarchies. It was also our intention to directly couple analytical with numerical 
solution methods to obtain an efficient alternative to other more standard formulations. 

2. Applied Goal: To use this new methodology to develop a system of analytical and numerical ejector/mixer 
nozzle models which require minimal empirical input. It is the goal of this research to provide more 
information than a control volume model; at a computational cost which is competitive with control 
volume methods. The DREA model is of direct use to the High Speed Civil Transport Program and is in 
the process of being adopted by both NASA and industry. Though most testing was done using 
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ejector/mixing data from the literature, new data was obtained using a hydraulic/gas flow analog model. 
These measurements provided information about the inviscid mixing interface structure. 


5.2.1. 1 General goal 

As described in Section 1, despite the rapid advances in both scalar and parallel computational tools, the 
large number and breadth of variables involved in both design and inverse problems make the use of 
sophisticated and even relatively simple (parabolized or boundary layer) fluid flow models impractical. With 
this restriction, one may concluded that an important family of methods for mathematical/computational 
development is reduced or approximate models. This model was related to other approximation methods as 
shown in Table 5.8 repeated here for convenience: 


Perturbation Method/Approximation 

Comments: 

Integral/Control volume 
(.i.e. 1-d or quasi- 1-d) 

Satisfy integral forms of conservation equations. Cannot 
predict local values. 

Current DREA model 
0(1) equation system 

Parabolic in terms of conservative flux quantities. Primitive 
quantities: u, v, p, p, M, locally predicted. PNS subsonic 
pressure approximation 

Boundary layer 

Parabolic in terms of u, v for all flows; p(x), thin layer, p(x) 
must be specified by an auxiliary relationship 

Viscous-shock layer 

Parabolic for supersonic flow; u, v, p(x) specified by auxiliary 
relationship, p(x,y) specified by inviscid y-momentum; PNS 
subsonic pressure approximation. 

Parabolized Navier-Stokes 

Parabolic (in space) for supersonic flow, elliptic subsonic 
regions; u, v, p(x), p(x,y) Viscous y-momentum. PNS 
subsonic pressure approximation 

Full Navier-Stokes Analysis 

Parabolic-hyperbolic, fully viscous equations. Always solved 
using unsteady or iterative method. 


TABLE 5.8 Classical and current approximation methods in fluid dynamics. 


To overcome this limitation, a successful combination of perturbation/numerical modeling methodologies, 
which provide a rigorously derived hierarchy of solutions, has been developed. Additionally, the solution to 
these models utilized analytical solutions to resolve singular behavior as required. Hence, classical methods 
are to be combined with efficient numerical methods to yield an efficient and original class of fluid flow 
models. 

By way of conclusion, it is the author’s opinion, that classical, asymptotic and reduced modeling methods 
are to be preferred over gross application of computational tools. Having said this, though, it is important to 
recall that perturbation methods and their intimate connection with classical methods have significant 
limitations, principally that they are difficult, if not impossible, to apply for problems where some form of 
linearization is unavailable, or where the nonlinearity is intimately connected to the problem. Progress in 
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handling highly nonlinear problems is discussed by Malmuth (1993) in a paper appropriately titled: “Some 
Applications of Combined Asymptotic and Numerics in Fluid Mechanics and Aerodynamics”. An interesting 
set of asymptotic analyses is provided in Ockenden and Ockenden (1995). 


5.2. 1.2 Applied goal 

In particular, the main objective of this research was to develop a system of analytical and numerical 
ejector/mixer nozzle models that require minimal empirical input. A code, DREA Differential Reduced 
Ejector/mixer Analysis has been written with the ability to run sufficiently fast such that it may be used either 
as a subroutine or called by an design optimization routine. We believe that this experimental validation of 
these models is provided by comparison to results obtained from the literature, industry data, as well as, 
dedicated experiment performed at Texas A&M. The error associated with the comparison of these models to 
the DREA simulation is summarized in Table 5.9. 


Experiment 

Mixer/ splitter 
plate 

relative error/comments 

Gilbert and Hill 
(1973) 

simple, 2-d 
slot mixer 

entrainment error 1.86% 

flow field prediction good, max. err. 8% 

Goebel and Dutton 
(1991) 

simple, 2-d 
slot mixer 

Flow field profile results: spreading layer 
width good, location, moderate 

Fernando and 
Menon (1993) 

Re VOIt =0.087 

h<A=0.428 

ho/h(0)=0.006 

1 

Flow field profile results: spreading layer 
width good, location, moderate. Failure to 
predict local acceleration at chutes. 

Amey and Lidstone 

simple, 2-d 
slot mixer 

Flow field profile results: spreading layer 
width good-moderate, location, good- 
moderate. CFG, 1 1%, entrainment 23% 
u design” NPR 

Thayer et al. 

Revo.,-2.0 

ho/h(0)=1.0 

CFG, 10% full range; entrainment 20-25% 
at “design” NPR. Failure to predict local 
acceleration at chutes. Reasonable 
comparison with “local” acceleration 
modification 


TABLE 5.9 Summary of results of comparison with test cases. 


It is possible to conclude that although the DREA simulation has limitations (especially in predicting 
entrainment for short shroud ejectors and prediction of local acceleration effects due to vortical mixing), it 
clearly provides a useable engineering predictive tool. Additionally, the DREA code uses a very limited 
number of empirical constants. In fact depending upon how one counts an empirical constant, there are really 
only two truly empirical constants, i.e. the effective viscosity and the local vortex strength used in the vortical 
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mixing analysis. On the other hand, the author must admit that at least some of the analyses developed here are 
so approximate that they are little better than order of magnitude or scaling arguments. Nevertheless, they do 
not require additional empirical input. Given this limited dependence upon experimental input, it is anticipated 
that the solutions will be valid (at least to the level of accuracy described here) for a wide range of problems. 
This discussion brings us back to the comparison of the DREA code with other models. Recall that the other 
models of at least computational speed approximately equivalent with the DREA analysis were, a code 
developed at MIT, Fung (1995), as well a code developed at Boeing by Clark (1995). The comparison table 
used in section 1, here Table 5.10 is presented: 


Criterion 

DREA (Texas A&M) 

Fung (MIT) 

Clark (CFA) 

Vortical Mixing 

yes; analytically based 
turbulence model 

yes, piecewise model 
and scaling functions 

no, multiple lobe 
only 

governing equations 

simplified 2-d, 
conservative 

scaling law 

l-d 

primitive variable 
profiles 

yes 

no 

two stream model 

empiricism 

limited 

moderate (“scaling” 
constants) 

limited (unknown) 

compressible 

expressly developed for 
compressible flow 

1-d and extensions 

l-d and extensions 

internal shock /choke 
structure 

yes; Fabri-choke and back 
pressure limited 

no, pumping only 

current version no 

computational speed 
(suitable for 
preliminary 
design/optimization) 

fast (user may choose 
control volume and/or 
mixing) 

fast 

fast 

experimental validation 

yes, literature, proprietary 
data and experiments 

yes, literature and 
experiments 

yes, experiments 

publications 

journal articles, 
conference papers, 
dissertation and reports 

thesis, conference 
papers 

conference papers 
(main code 
proprietary) 

status availability 

preliminary summer 
1996, fully tested and 
documented 05-97 

complete 

complete 


TABLE 5.10 Comparison of mixer/ejector models suitable for preliminary design/optimization. 


One may reasonably conclude, that the DREA code does indeed provide a significant improvement over 
the other current models all of which have been developed to fill need for a fast, yet sufficiently powerful 
computational model. 

Finally, then an analysis code, DREA, i.e. one for which geometry and operating conditions have been 
specified has been developed. Performance characteristics (thrust, noise, degree of mixing, and pumping) are 
predicted in an efficient manner by this code to an adequate level of accuracy. Ultimately, this type of tool may 
be used to generate off line performance information in the form of "ejector" maps, which are tabular lists of 
performance characteristics. It is also sufficiently computationally efficient to be run in a "real time fashion or 
serve as the basis for an inverse design type of tool. 


5.2.2 Recommendations for further work 
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In this section modification to the current work are discussed. As with the conclusion, we split our 
recommendations into two types: 

1. Near term modifications to the DREA, 0(1) system of equations. These types of modifications 
are intended to leave the basic solver and governing system intact, while improving the 
comparison with experimental data or extending the range of validity of the current tool. 

2. Long term, major modifications to DREA code. These include modifying either/or the basic 
solution method or the fundamental formulation. 

Of course, these two areas of recommended further work are not mutually exclusive. For example if one were 
to upgrade the mixing equations to model wave effects, it would make little sense not to simultaneous upgrade 
the entrainment model, which ultimately provides the initial conditions to the marching equations. 


5.2.2. 1 . Recommendations for near term further work 

In this section, a brief discussion of potential moderate extensions to the current research is presented. 

These include: 

• Improvement of ideal and non-ideal entrainment computations for both back pressure dominated and back 
pressure independent (aerodynamical ly choked) flows. Though it is unquestionably possible to 
incrementally improve and upgrade our simple ordinary differential equation models, which estimate 
losses in entrainment, or even invoke an iterative matching procedure using our full equation, the fact 
remains that modeling flows with large regions of subsonic flow will be difficult using single pass space 
marching, i.e. “PNS like” methods. Examples of flows like this are presented in Figure 5.31. 

• These flows are fundamentally elliptic, demanding an iterative or relaxation solution technique. This area 
remains a challenge for CFD development work. A potential opportunity may exist for our combined 
analytical/numerical methodology, where by, regions of the flow field, such as the elliptic portions may be 
modeled using analytical tools, while well posed supersonic regions may be modeled (in their full non- 
linearity) using numerical methods. 
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FIGURE 5.31 A schematic of flow fields with substantial subsonic flow regions. 


• Developing a simplified “map” of basic ejector characteristics. This is a difficult process, in that there still 
does exist a completely adequate physical criterion for predicting the ejector operation mode (i.e. subsonic 
exit or supersonic flow at the exit). Analytical and control volume methods give us perhaps the best 
chance of developing a rational basis to understand the scope and limits of ejector operation. 

• Develop a methodology to model truly three-dimensional vortical effects with local acceleration effects. If 
one can accept a moderate development effort, these effects may be modeled using a new family of 
approximations. Consider the three dimensional mixing model: 


cty 

3c 


V eff 

U 


d <f> 

& 


eff 


UJ 


S' <f> 

dz 1 


(5.5) 


The initial conditions associated with this type of flow are shown in Figure 5.32. 
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FIGURE 5.32 Geometry definitions for a single lobe, 3-d model. 


Unfortunately, there is no simple decomposition method that will permit us to model this three- 
dimensional flow in a completely satisfactory manner. To solve a relationship, like the one found 
above, it is necessary to develop a three dimensional solution, i.e. <(>=<|>(x,y,z). 

Though a solution is not developed in significant detail (since it outside the scope of the current 
research), a purely analytical solution based upon 2-d Fourier transforms (alternatively Green’s 
function expansions) is presented. The governing equation solved (note the local “y” and “z” 
linearizations): 


d<f> 

3c 


, d <f) 

= Cl xyX ~y 

dy 


+ a xzX 


d l (f> 

dz 1 


(5.6) 


with the boundary conditions described in Figure 5.32. The solution is: 
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(5.7) 


Equation (5.7) is a double infinite series in terms of Green’s functions. Like it’s 1-d counterpart, it 
seems reasonable to expect, that this solution would be very efficient, i.e. requiring a minimum number 
of terms, close to the initial condition. Extension of the numerical solution to 2-d is significantly more 
difficult. Here it would also be necessary to develop some form of a higher order solution method 
appropriate for parabolic equations, for example a fourth order, compact ADI scheme, see Hirsch 
(1975). 
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• Continuation of testing of the current DREA model exploring the possibility of a static pressure correction 
equation. This would employ a decomposition of the form: 

p = COp + (1 - 0))p (5.8) 


where the term co is the Vigneron (Vigneron et al. (1978)), (See also Appendix E)) parameter 
developed earlier and used to provide a well posed spatial marching method for subsonic regions. Its 
value is: 


yM 2 

(0= T 

1 + (/ - \)M 2 


(5.9) 


for M<1 and 1 for M>1. Using this parameter, the pressure is split into two portions: (1) a locally 
computed value and (2) a boundary layer portion (a function of x only). Of course, this decomposition 
requires a knowledge of the Mach number distribution, M(x,y), which is not explicitly known. As such 
it is not a completely satisfactory solution. 

• Consider the introduction of inverse design methodologies. Clearly the reason to develop any a design 
code is to have a tool of sufficient efficiency such that trade studies are possible. Automating these studies 
using optimization methods (see Luenberger (1984) for classical gradient based methods, and Goldberg 
(1989) for genetic based methods or simulated annealing (Press et al. (1992)) is course, desirable. The 
computational cost of identifying global minima can be truly amazing. Upon considering the cost of 
performing a direct search in say N variable (of course this is highly inefficient, but it can and will 
delineate a local minima) using a 100 node mesh, one is now faced with performing at least 100 N marching 
computations. For N=3 and assuming 6 seconds/ run it would be expected to require a PC for on the order 
of 70 days! Of course, we would probably not use such a conservative algorithm, but optimization is a 
dangerous (and expensive) game. Additionally, only running the code for a control volume case would 
reduce the cost by half. 

This discouraging assessment of optimization costs is offset when one considers the use faster scalar 
machines or indeed parallel machines (See Geist et al. (1996)). A message passing variant of this code 
and an optimization algorithm (for example a genetic method) could make robust optimization extremely 
competitive. Of course, this places a tremendous impetuous upon us to develop the most realistic, accurate 
and reliable basic model possible. 


5.2.2.2. Recommendations for long term further work 

In the realm of long term recommendations, the possibility of considering new fundamental governing 
equations and modeling effect heretofore unattainable using more limited equations is considered. Clearly, the 
most pressing modeling concerns for the governing equations system is: 

1 . Directly coupling wave effects 

2. Dealing with the static pressure term in as a complete manner as possible. 

3. Modeling three dimensional effects. 

We address these limitations next. 

Section 1 developed equations that were even less complex than the boundary layer family . Unfortunately, 
it was shown that, though the mixing was well modeled, the wave phenomena, such as, expansion 
compression/shock systems and normal shocks could not be directly modeled at all. To deal with these 


171 



features, a double valued transformation was used, which contained two possible solutions. This worked well 
in the streamwise direction, but lead to significant problems in the cross-stream. To deal with this situation, the 
modified inversion transformation was introduced and used in transonic mixing regions (recall this required 
neglecting a portion of the static pressure field). Additionally information describing the location of the choke 
was also needed. These limitations clearly indicate, in spite of the fact that the equations used were very simple 
and linear, cannot serve as a completely satisfactory model. 

With this in mind, it is of interest to consider a set of equations that explicitly model two-dimensional 
wave phenomena and mixing simultaneously. This set is rather easily written: 

x-momentum: 


d , d dp d 

— (Pu) + — (pv)+- — 

dx dy dx dy 


d_ 

dy 


di k 


(pU'^i-iPeff^) 


y-momentum: 


4- (pw) + 4~(pv 2 ) + ^r = 0 

ox dy dy 


energy: 


3 / t/1 , 9 / u\ 9 / >U>\ 9 

— (puH) + — (pvH) = - — (pv'H') = — (- — ) 

dx dy dy dy Pr, urA dy 


and continuity: 


4 ~( p u ) + 4 -( p v ) = o 

dx dy 


(5.10) 


(5.11) 


(5.12) 


(5.13) 


Note that as before, that the Prandtl number has been assumed to be close to unity, hence dissipation terms in 
the energy equation, may be neglected 

Given suitable turbulence modeling for the right hand sides (mixing terms), it is claimed that this set will 
support both viscous mixing (parabolic) and hyperbolic wave expansion. This claim is verified by the fact, that 
given isentropic flow (irrotational), and neglecting the viscous and heat transfer terms, that we derive the 
classical (non-linear) potential equation: 


[i-o;]4>„+[i — ^=o (5.i4) 

where O, is the velocity potential and a= speed of sound. The speed of sound is related to the energy equation, 
by noting that H 0 =constant satisfies energy equation and using state. Unfortunately this equation set equations 
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(5.10)-(5.13), which are essentially a set of “viscous shock layer equations” are computationally difficult to 
deal with. 


The question arises: under what circumstances would the perturbation expansion system developed in 
Section 1 first exhibit explicit wave phenomenon? To ascertain this one must consider the absolute minimum 
governing equation system required to derive a linear wave equation. This minimum system is written: 

x-momentum: 


pu 


du dp 
— + — 
3c 3c 


= 0 


(5.15) 


y-momentum: 





continuity: 

d{pu) | o{pv) _ Q 
3c dy 


(5.16) 


(5.17) 


energy: 


ojpuH) | c{pvH) _ q 
3c dy 


(5.18) 


a form of the state equation: 


dp = -a 2 dp 


(5.19) 


and finally, the irrotational condition: 


3c _ du 
3c dy 


(5.20) 


Clearly, it is apparent that given the absence of mixing terms, that the governing equation set is complete even 
when considering only the 0(1) mixing problem, except for the final irrotational condition. It is this condition 
which, according to our expansion suppresses wave phenomenon. Considering this last equation, we note that 
even with “y” scaling at 0(1), that the first term in our expansion for v, i.e. v 0 is assumed to be O(e). 
Therefore, wave compression/expansion features will be suppressed till O(e). To emphasize this, we write the 
0(e) wave equation system: 
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= other terms 


(5.21) 



jl ^ ; o 

1 3c dy 


and 


3> o du 2 

= 1 - other terms 

3c dy 


(5.22) 


where the “other terms” signify other expansion terms which do not contribute to the structure of the linearized 
wave portion of the equation. 

Hence, it is apparent that wave phenomenon have been suppressed to higher order using the perturbation 
expansion developed here. The order of magnitude limitation (and the one common to all boundary layer 
analysis or parallel flow analysis) is associated the cross-stream velocity, v, to be 0(e). This prompted the 
supplemental analyses, i.e. the inviscid choke location described in Section 2, and the discussion of more 
complete (for example the viscous shock layer equations) developed previously. 

Finally it is necessary to consider 3-d effects. If we are willing to accept a family of governing equations 
of considerably greater complexity, such as the VSL equations, modeling of three-dimensional effects becomes 
relatively straight forward. Our governing equations take the form: 

x-momentum: 


d_ 

dx 


d_ 

dy 


d_ 

dz 


dp 

dx 


d_ 

dy' 


du, 


(pu’) + ±(pm’) + ± (pun) + -y - = | (/*„ 


du. 


(5.23) 


y-momentum: 


d . . d . 7 . d . dp 

— (puv) + — (pv 2 ) +—(pvw) + — = 0 (5.24) 

ox dy dz dy 


z-momentum: 


d . d . d . 2 1 dp d dw d dw 

— (puw) + — (pm) + — (pw ) + — = — (Peff —)+ — (Meff — ) 
dx dy dz dzdy dydz dz 


(5.25) 


energy: 


±(puH)+^-(pvH) + 

dx dy 


d_ 

dz 


(p»H)-4-(~y 


dy Pr, 


iurb 


3Hf, | d 

dy dz ?r lurb dz 


(5.26) 


and continuity: 
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(5.27) 


d . . d . . d _ 

—(pu)+—(pv) +—(M = Q 

ox oy 02 

Once again, the Prandtl number has been assumed to be close to unity, hence dissipation terms in the energy 
equation are neglected. 

The expansions are more complex, since the unknown “z” direction mass flux: pw, has been added to the 
solution set and the resulting equation system will need to be factorized as described in Anderson et al. (1984). 
It is important to point out that in development work it is necessary to be cautious not to “reinvent the wheel” 
since the equations that we have now described are indeed essentially PNS equations. With this complex set of 
equations, we have come almost full circle. As such the development of modeling methodologies described 
here is concluded. 
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APPENDIX A 


FORMAL HIGHER ORDER EXPANSION, 0(e ,/J ), EQUATION SYSTEMS 

As discussed in Chapter I, one of the main strengths of the overall modeling methodology is development 
of a formal perturbation hierarchy. This type of expansion method permitted the "apriori" estimation of the 
error associated with a particular approximation, as well as, the formulation of corrections to the lower order 
approximations. In this appendix, formal development of both the first order, i.e. 0(1), and second order, 
0(e i/2 ) expansions are derived. Since the first order expansion has been developed in the text, concentration is 
placed upon the mathematical character of the second order expansion. 

Our work in this section will be to develop a system of partial differential equations describing viscous 
mixing. As such, we begin by considering compressible, two-dimensional Reynolds equations. The Prandtl 
number has been assumed to be close to unity, hence dissipation terms in the energy equation, may be 
neglected. As before, the conservation equations are written: 

x-momentum: 

4~(Pu 2 ) + -^-(Pv 2 )+^- = -~-(pu'v')-^-(pu'u') (A.l) 

dx dy ox ay ox 


y-momentum: 

(puv) + ^(pv 2 )+^ = --^(pu' v ')--^-(P v ' v ') ( A - 2 ) 

dx dy dy dx dy 


energy: 

< A3 > 


and continuity: 

(pu) + ^-(pv) = Q ( A - 4 ) 

dx dy 

To proceed with the analysis, it is first recognized that when non-dimensionalized, there are a number of 
"small terms” inherent to the governing equations. The small terms come from both the magnitude of the 
dependent and independent variables. For many problems, the flux difference between streams is relatively 
small, as is the mixing layer thickness compared to the streamwise distance. Formally stated this yields a set of 
characteristically small parameters: 
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where 6^ refers to the mixing layer thickness. As indicated in Chapter I, it should be noted that e«e G may not 
be small for all problems of interest. It will be shown, though, that this choice of expansion parameters 
provides reasonable results despite the fact they are formally not good approximations. Since it is important to 
have differential equations that can resolve the mixing layer itself, it is necessary to rescale (stretch) the cross- 
stream or "y M coordinate: 


L * 

x * 0 ( 1 ) y = £ 2 y 


d _ 1 d 

*y m ei*y 


As before, the streamwise variables are simply expanded: 


p/ = p/o + Pf]s 2 + pf' 2 £+... f = u\ H 
g = So + Si £~2 + g 2 e+... g =p , u , H* 


(A.6) 


(A. 7) 


Recall that the expansion of the cross-steam flux and velocity took more care. It is expected that for the 
basically parallel flows of interest, that the cross stream velocity, v, is very small, 0(e). The cross stream mass 
flux, pv, is also small, but in order to conserve mass this flux must be an order of magnitude larger, i.e. 0(e I/2 ). 
With these estimates, the cross-stream mass flux and velocity are written: 

* • ^ . * a. 

Pv - PV 0 S 2 + 

(A.8) 

3 

V*~V o£ + V/£ 2 +V2£ 2+ **- 


Proceeding with the scaled problem, the magnitude of the Reynolds stress terms (note we have suppressed the 
* denoting non-dimensionalization of the dependent variables to simplify the notation) are estimated. 
Following the development given in chapter I we give to 0(e I/2 ) the conservation relationships: 


pu f v' ~ - 


U , 


G s 2 


\*y 


(pu 0 +O0W,H o + pu 0 U l )s 2 +...\ 


(A9) 
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dy 


~.(pu Q )e 2 +. 
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Notice the last term which resulted from our linearization of pG’u’/u. This non-linear term will provide the 
fundamental coupling between the 0(1) system and the 0(e l/2 ) system. 

Since the interest in this appendix is in higher order expansions, it is necessary to estimate the size of the 
other Reynolds stress terms. Consider the other Reynolds stress term pv’v’ (which is the transport of turbulent 
cross-stream momentum flux by cross-stream turbulent velocity): 


pv V 


JL 
U . 




-(puv) = - 


U) 


G e 2 


\*y 


u oP v o £ 


(A. 10) 


The term, pu’u’ (which is the transport of turbulent streamwise momentum flux by streamwise turbulent 
velocities) is estimated to be: 


pu'u'' 


l c/I 


£ g-t(p uv ) = - 

ox 


U 


V eff 1 ^ 


£G fa \. PU ° V ° £+ - ”] 


(A.ll) 


which is 0(s 2 ) and, as expected, is much smaller, than either pu’v’ or pv’v’. 

Using equation (A.9) and invoking the additional well founded approximation that the turbulent Prandtl 
number Pr t =0(l), the Reynolds transport of total enthalpy is written: 


pH’v’ = 


1 

(v A 
eff 

1 

o 

Vi 

« 

Pr, 

l u ) 



pu 0 H 0 + (pu^Hq +pu 0 H i )s 


(A. 12) 
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dy 


Substituting the expansions and collect terms yields, the two lowest order systems. As before, the 0(1) 
system may be written: 

x-momentum: 


— (P“l + Po) = 
ox 


d 

dy 



G 


£g d 
e dy * 


(P“l)\ 


(A. 13) 


y-momentum: 


187 



(A. 14) 


d£o 

dy 


= 0 


energy: 


OX 



Pr. 



£ dy 


(pu 0 HJ] 


and mass conservation: 


(A. 15) 


d , . 

— (pUo) + 
ox 


~(pvj = 0 
dy 


(A. 16) 


This is the 0(1) system, that is solved. Again, as it is written, it is not possible proceed, since the last term in 
the preceding equations is not closed. To solve this system it is necessary to approximate the cross-stream 
mass flux term: 


]_ 

pv 0 £ 2 


Sc, 


V _JJ_ 

u , 


dy 


~.(pu 0 ) 


(A. 17) 


where Sc t is the turbulent Schmidt number, a 0(1) term. Equation (A. 17) is essentially a linearization of our 
Reynolds stress closure. 

One more convenient modification can be made. Since the cross-stream pressure gradient 5p/5y«0(e) 1/2 , it 
is permissible to add this term to the right hand side of the streamwise momentum equation thus yielding: 


d , 2 

— (pu 0 +Po) = 
ox 
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£ g d 
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(A. 18) 


This completes the formal perturbation development of the system 0(1) of equations used in this project. 

Now, the 0(e 1/2 ) system is also presented in this appendix. The conservation equations take the form: 
x-mo men turn: 
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y-momentum: 


and continuity: 
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energy: 


d 3 

— (pu 0 H \ +P“iH 0 ) + — (pv 0 H 0 ) = 
dx dy 
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dy Pr, \U ) Q s dy 
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+ 0 <A ' 22) 

dx dy 

A state equation, which is assumed to be a thermally perfect ideal gas, is included. Combining the ideal gas 
state equation with the definition of total enthalpy and retaining second order terms yields: 


P\ = 


i +{p“o u i + 

7 



(A.23) 


Equations (A.19)-(A.23) provide the governing set of equations to 0(e 1/2 ) which are listed in Table A.l. 
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Equation 

Unknown 

momentum, equation (A. 19) 

U| 

energy; equation (A.21) 

H, or Ti 

mass; equation (A.22) 

pi 

y-momentum equation (A.20) 

pvi 

state; equation (A.23) 

Pi 

Total=5 

Total=5 


TABLE A. 1 Equation versus unknown count for 0(e 1/2 ), system. 


Comparison is now made between the 0(s ,/2 ) system with the 0(1) set and other classical fluid flow 
approximations. The 0(e 1/2 ) family of equations is essentially linearized boundary layer or, alternatively stated, 
the 0(e 1/2 ) system are essentially linearized second order boundary layer equations (Van Dyke (1975)). The 
most striking addition to the 0(e l/2 ), is that now, the y-momentum equation, equation (A.20) has become non- 
trivial. As such, pi or pv t does not require the addition of the mixing mass flux closure relationship. Similar to 
classical boundary layer approximations, though, the cross-stream pressure gradient is explicit. In our system 
P! is computed using explicit information from the 0(1) solution, through equation (A.20). This bears close 
analogy with classical second order boundary layer equations where the pressure gradient is coupled to the 
streamline displacement using the inviscid solution of Lighthill (1958). Displacement thickness and the 
pressure gradient may be directly computed in our 0(e i/2 ) problem by invoking local continuity and local cross- 
stream momentum conservation. As before, other quantities, such as Mach number and total pressure, are 
supplemented by their own specific definitions. 

This completes the formal perturbation development of the systems of equations used in this project. 
Recall that the goal was to obtain governing equations that were simpler than classical boundary layer, but 
more powerful than quasi- 1-d equations. This goal has certainly been achieved by the 0(1) system. 
Additionally, the higher order 0(e 1/2 ) system has been shown to retain many of the characteristic of classical, 
second order boundary layer methods. This close connection is to be expected, since the higher order 
equations should be consistent with the classical perturbation hierarchy, Table 1.3; a portion of which is 
repeated below as Table A.2. 


Perturbation Method/Approximation 

Comments: 

Integral/Control volume 
(.i.e. 1-d or quasi- 1-d) 

Satisfy integral forms of conservation equations. Cannot 
predict local values. 

Current DREA Model 
0(1) equation system 

Parabolic in terms of conservative flux quantities, primitive 
quantities: u, v, p, p, M, locally predicted. PNS subsonic 
pressure approximation 

Boundary layer 

Parabolic for all flows; u, v, p(x), thin layer, p(x) specified 
by auxiliary relationship 


TABLE A.2 Classical perturbation approximations and current DREA model. 
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This concludes our discussion of the higher order perturbation system. It is worth noting however, that 
both the previous systems are essentially parabolic in nature. As discussed in Chapter V, the wave effects (i.e. 
expansion-compression) structures inherent to our problem are not be modeled by either the 0(1) nor the 
0(e 1/2 ) equation sets. To do this requires using our expansions, retaining terms through 0(e 3/2 ). The 
considerable delay in the appearance of wave effects indicates, that to model wave phenomena, that it may be 
necessary to start with an entirely different perturbation expansion scaling. 
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APPENDIX B 


CRITICAL BACK PRESSURE COMPUTATION AND SUBSONIC FLOW LIMITING 

CONDITIONS 


The text developed in some detail two fundamental modes of ejector flow, namely, back pressure 
constrained and back pressure unconstrained flows. These two flow regimes have significantly different 
modeling methods associated with them. Although the distinction that two flow problems exist is clear, at no 
point did we discuss how one chooses, a priori, which is the appropriate model. In this section an approximate 
model to identify the appropriate flow regime is discussed. In addition to the supersonic/subsonic regime 
transition there are a series of constraints imposed upon the subsonic operation of an ejector. Attempting to 
operate outside of these constraints will cause DREA to fail due to physical reasons. As such, it is of 
considerable interest to provide an analysis to delineate these operational regions. Though it is beyond the 
scope of this project to develop operational maps using the full suite of constraints, as an example, one of the 
“back flow” constraint problem is modeled. Results of the back flow problem are compared with DREA to 
show that the convergence failure of DREA is indeed associated with one of the physical constraints, rather 
than some numerical difficulty. 


B. 1 Critical back pressure computation 


The choice of the proper condition is dependent, as one would expect, upon the back pressure or exit 
pressure level that the ejector senses. Here, a value for the exit pressure which separates the back pressure 
dependent regime from the back pressure independent regime is developed. This value is defined as the critical 
back pressure or p crit . The computation of this value follows. Referring to figure B.l, we consider the 
combined model. 


M=1 

Aerodynamic, "Fabri" choke 



P=P 

crit 


M 


3 


FIGURE B.l Control volume definitions for the critical back pressure computation. 


The first portion of the model involves the Fabri limited flow. The solution for this problem (which is 
described in Section 2.2.4 are the conditions associated with supersonic flow and choking in the secondary 
stream. Hence, if the exit conditions associated with this flow are known, then so would be the critical back 
pressure level. 

Given a value for the critical back pressure, the two possible flow regimes are defined: if the back pressure 
level is below p crit , then the Fabri-choke is the correct solution. Alternatively, if the back pressure level is 
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greater than p crit then the flow is not properly modeled using Fabri theory and should be modeled using the 
back pressure dependent model. These criteria are summarized in table B.l 


( Pe )< ( Pcnl ) 
Poi Po^ 

Back pressure independent flow 
Section 2.2.4 

‘I 3 

Al 

«*| s 

Back pressure dependent flow 
Section 2.2.3 


TABLE B.l Critical back pressure formulations. 


From Table, B.l it is clear that the key parameter governing the ejector flow mode is the value of p C ri/Poi- 
It is possible to compute p crit /p 0 i by reformulating a set of 1-d conservation relationships. These are repeated 
for convenience: 


mis + rii2s ~ m 3 


(B.l) 


misUls + P Is A m2sU2s + P 2s A2s - m3U3 + P 3 A3 


(B.2) 


misT ois + rii2s T 02s - m 3 To3 


(B.3) 


and the exit pressure definition: 


P3 = Pc =P«, 


(B.4) 


Combining equations (B.1)-(B.4), it is possible to obtain a relationship in terms of M 1$ , M 2s , and the critical 
pressure ratio p cr i/poi- Provided below are the details of this procedure. Starting by rewriting the secondary or 
primary stream inlet terms as functions of M 2s or M| S yields: 


m = puA = 



r-i 

2 


M 1 


r_± i 

2(1— r) 


(B.5) 


mT = 


yMP 0 Aj 


0 r„ 



Lzl 

2 


M 2 


r± i 
2(i-r) 


(B.6) 
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mu + pA = P 0 A 


(B.7) 


1 + 


r-i 

2 


M 2 




To proceed, it is necessary to rewrite the exit conditions in terms of the upstream quantities and the exit 
pressure p crit . Mass and energy can be easily rewritten in terms of upstream quantities, while, momentum is 
more difficult due to the velocity term (a non-conservative variable). To eliminate the exit velocity, it is 
possible to use energy and state: 


m/s Hois + ri/2s Ho2 S = E = 



Pcr,,A U 3 

w 



(B.8) 


where W=m 1 +m 2 . This is a second-degree equation for u 3 . Solving it yields: 


l 



( r Peri, A 1 

2 1 

I 2 E 

j 

7 Peri, A 


I7-1 w ) 

1 w 

1 

1 


(B.9) 


Recall that E and W are known in terms of the upstream variables. Rewriting momentum: 


misUl, + Pi s Als + m 2 sU 2 s + P 2s A2s = m3U3 + Peri, A 3 (B.10) 


where (using energy) mu 3 is written: 


]_ 

2 

+ 

(B.ll) 

Per, A 3 

r - 1 


mu y = 


r 


r-i 


Per, A 


+ 


2yR 


r 


~{mis Tois + m 2s To 2 s){m u + m 2s "j 


Now, normalizing equation (B. 11) by PoiA h and substituting our expressions for the upstream quantities 
in terms of M 2s or M ls yields (with considerable manipulation) the non-linear, algebraic equation: 
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(1 + — Mu)h (1 + r M 2 ,J + (—)(—)(! + Ml') >-r(l + Y Ml) - 

2 p 0l A, 2 

+[— (~) 2 + -[Ml(l + Zy- M 2 ,s0 + 

Y-i Pot r-i 2 

+ Mu Mu (^)A) <1 + ^->(1 + tj- M’„0r 1 a + ^ A A>wr, + 

p 02 A, T 02 T 02 2 2 


(B.12) 


+mU— /(— Y d* 

Pot A, 


Y- } 

2 


y+J l 

Mis) i-r f 2 J 


+ 


_L(Pj-)(Al ). 0 

Y-l Poi Ai 


This rather cumbersome expression, equation (B.12) relates the critical pressure ratio p C n/poi to the 
aerodynamic choke quantities, Ai s , Mi s , A 2 S and M 2 S . Solution of these relationships is relatively straight 
forward, in that, p cn /Poi is decoupled from the Fabri analyses, permitting solution of the Fabri problem and 
inversion (using a single variable, nonlinear equation solver) to compute the critical pressure ratio. 


B.2 Subsonic limiting flows 

In this section the limiting flow field conditions for subsonic mixing are discussed. The meaning of a 
“limiting” subsonic flow in this application describes flow fields for which the intended ejector operation, 

i.e. primary stream flow inducing secondary entrainment, thereby, satisfying the downstream pressure 
constraint fails. It is possible to identify three possible limiting subsonic limiting flow situations: 

1. Incipient reverse flow into the secondary inlet. 

2. Choked flow in the secondary inlet. 

3. Choked flow in the exit mixing stream. 

These three cases are illustrated in Figure B.2: 


Potential backflow condition or 
secondary stream choke 



Potential mixed stream choke 
FIGURE B.2 A diagram showing possible limiting subsonic flow fields. 
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Our task is to determine upper and lower bounds on (P 02 /P 01 XA 2 /A 1 ), such that there is sufficient secondary 
pressure to prevent back flow, without inducing secondary inlet choking or exit stream choking. Hence a 
(P 02 ^P o i X A 2 / A j ) cr it used to delimit various flow regimes is to be developed. 

Analysis of the critical condition problems may be performed by considering the system: 


rh x + m 2 = m 


(B.13) 


mi ui + P/As + rii 2 u 2 +p 2 A2= riii ui + p 3 Ai 


(B.14) 


riiiToi + rii 2 T 02 — riii T 03 


(B- 1-5) 


and the exit pressure definition: 


P 3 = P« 


(B-16) 


Equation (B.13-B.16) can be placed into a more useful form for characterizing the critical flow conditions 
by rewriting them: 



(B.17) 


(i -V+— 
2 2 


M\ +[(2-G)/]M 3 2 +1 = 0 


(B- 1 8) 
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Ml 

( p \ 
r 0\ 

(a] 


1+ f*f|(V 
v JWj J V ^oi > 

1 ■ 

2 

h • ^ 

i + 

V m x ) 

1 

2 

mJ 

VpJ 

uj 

M 

0+V*) 

fl + 7 1 Aff 

V 2 

r+i 

J 2(r-i) 


(B.19) 


Completing this set of equations are the definitions: 





(P 1 
■*02 

hi 

f f \ 
J 02 

\M X ) 

UJ 


UJ 


1 


+ 


r-i 

2 



K 


r_± i 

|2(r-D 



(B.20) 


Hi 

«i 


f M 0 

(t \ 
*02 

l MJ 

k T 0l J 


1 + — — - A/, 2 

^ * 


1 + 




]_ 

2 


(B.21) 


Equations (B.17-B.21) provide the governing equations for the mixing problem in terms of the quantities of 
interest. As such, we may now pose the three limiting flow problems, listed in Table B.2. 


Flow Limit 

Unknowns 

M 2 =0 

M„ (Poj/PcXAj/A.U, g 

m 2 =i 

M 3 , (Po2/Poi)(A2/A,)cri., G 

m 3 =i 

Mb (Po 2 /Poi)(A 2 /A,) cril . G 


TABLE B.2 Limiting flow problems. 


where the parameter of greatest interest is, again, (P 0 2 /PoiXA 2 /Ai) C nt. The above problems may be solved 
using any one of the nonlinear system solvers developed previous. 

An example of one of these limiting flows is an ejector which has a sufficiently high back pressure (one 
of several parameters which might cause this) such that the flow does not exit from the mixing duct exit, but 
reverses in the mixing duct and dumps out the secondary inlet (the first case in Table B.2). This flow field 
is illustrated in Figure B.3: 
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FIGURE B.3 A limiting flow problem, i.e. potential incipient back flow. 


This type of operation is clearly unacceptable, both physically and mathematically. Our task is to identify 
conditions where this flow is imminent and to show that this is where the DREA simulation code becomes 
unusable. 

The condition for incipient reversed flow into the secondary inlet is relatively simple; namely the 
secondary velocity and mass flow rate are zero, i.e. U2=m 2 =0. Hence, the governing equations reduce to 
relatively simple forms. The unknown in this system is ultimately (Po 2 /PoiXA 2 /Ai) crit . This parameter is 
important because if the ejector operating conditions have (Po 2 /PoiXA 2 /A 1 )>(Po 2 /PoiXA 2 /Ai) cr j l > then back flow 
will be prevented. Conversely for (Po 2 /Poi)(A 2 /A I )<(Po 2 /Poi)(A 2 /A 1 ) c rit, reverse flow will occur and the DREA 
solution method is ill-posed: 



REW 



(B.22) 



+ 


Gy 


Ml +[{2-G)y]Ml + 1 = 0 


(B.23) 


Ml 

Mj 


■ 01 


PJ 


f A \ 


UJ 




- 1=0 


(B.24) 


Equations (B.22-B.24) may be solved to yield the critical value; i.e. (Po 2 /Poi)(A 2 /A 1 ) crit . As an example 
of the use of this method; consider a subsonic constant area mixing problem. The application of DREA for a 
primary NPR value of approximately 2.09 ends in convergence failure. The DREA total pressure area ratio, 
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(P 02 /Po,)(A 2 /A|)=8.76. It is suspected that this convergence failure is due to incipient back flow into the 
secondary inlet. Is this true? Solving equations (B.22-B.24) at various primary NPR’s and comparing with 
the DREA simulation yields Figure B.4. 


DREA Code limiting run 



o 

Oh 


5 1 

1.8 1.9 2 2.1 2.2 2.3 2.4 

Primary Stream NPR 


FIGURE B.4 Comparison of a limiting flow solution with a non-convergent DREA simulation. 


From Figure B.4, it is apparent, that the convergence failure in the DREA simulation is indeed caused by 
incipient back flow into the secondary stream. 

Unfortunately, due to the relative complexity of the necessary compressible formulation, it is difficult to 
get a sense of the influence of the parameters. As an aid to understanding, it is useful to consider an 
incompressible solution to the critical pressure ratio problem (this is effectively the secondary pressure needed 
to prevent reversed flow): 


M 

hi 

vpj 

v A/ 


p« a l r a i 
Pi A Pi vA J 


(B.25) 


From equation (B.20) the following trends may be identified: 

1. Large back pressure; i.e. small NPR causes an increased likelihood of back flow. 

2. Large exit stream area causes an increased likelihood of back flow. 

3. Large primary mass flow (i.e. large u,) causes an increased likelihood of back flow. 

Though these conditions are much more easily ascertained using this incompressible analysis; it is believe 
that the trends identified using equation (B.25) will remain consistent even for the compressible problem. 

This completes our discussion of limiting subsonic flow cases as well as subsonic/supersonic transition 
flows. Ultimately an operational map of limiting criteria is to be designed but is beyond the scope of this 
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study. As a demonstration of the ultimate importance of estimating operational characteristics, consider a 
LERD restricted, subsonic ejector nozzle which is presented in Figure B.5: 



The upper and lower bound upon primary stream pressures are immediately obvious from Figure B.5. As a 
first step to obtaining the types of maps, it may be noted that the governing equations provided in Appendix 
D provide the preliminary basis for these operational maps. 
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APPENDIX C 


INVISCID SUPERSONIC/SUBSONIC INTERFACE ANALYSIS 

As discussed in Section 2, estimates of the boundaries of subsonic and supersonic flow are required to 
model the classical supersonic ejector problem. These transitions between flow regimes are directly related to 
expansion-compression effects in the ejector flow field. A simplified analysis is developed here to provide 
estimates of these regions. 


C. 7 Problem description 

The problem under consideration is a two-dimensional, unconfined, supersonic primary stream with a 
subsonic secondary stream shear layer as shown in Figure C.l . 


free stream 
(2); U , M<1 


Xc/2 first critical 
< > 


nozzle radius, y=l 


w 


u 


interface slipline 


(1); M , M >1 


FIGURE C.l Diagram of the unconfined supersonic-subsonic interface (slipline) problem. 


The following assumptions are made: 

(1) The flow quantities may be decoupled into a base quantity plus a perturbation value. Non-linear 
terms are neglected. 

(2) The fluids are ideal gases. 

(3) The flow is inviscid and irrotational except at the interface, where h=r=l. 

According to Pai (1954) this is an acceptable approximation for large Reynolds number flows at locations 
where turbulent mixing is still relatively unimportant. Further, the vortex sheet approximation is a fundamental 
assumption of the Fabri choke analysis (Fabri and Seistrunck (1958) and Addy et al. (1981)). 

The flow is modeled as irrotational, so a small disturbance velocity potential, ((►', is used: 


Ui = Ui + Ui = Uj + 


dx 


(C.l) 
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where the subscript M i" describes the stream, primary (1) or secondary (2). For later reference, the governing, 
small disturbance, potential equations, (Anderson (1985)) for the primary (supersonic) stream are: 


„2 d% d% 

Pl d x 2 ' dy 2 


= 0 




and, for the secondary (subsonic) stream: 


Pi d x 2 dy 2 


= 0 


;fi 2 2 = l-M 2 2 


(C.2) 


(C.3) 


This completes our description of the governing equations. In the next section a brief review of the 
classical results giving several important relationships for later reference is presented. These results are limited 
in applicability to small stream pressure differences. Afterwards a description of the formulation of a strained 
coordinate modification to the slipline model, applicable for a larger range of static pressure ratios, is 
provided. 


C.2 Review of classical results 

This section discusses the strength and weaknesses of two classical, approximate methods which include 
an inlet perturbation model and a slipline perturbation model with its associated base flow model. 


C.2. 1 Inlet perturbation model 

This section describes an inlet perturbation model for confined flows which linearizes the flow about the 
inlet stream quantities. The analysis is closely related to normal mode methods and, as such, contains no initial 
condition information. Additionally, static pressure imbalances must be very small. From Pai (1952) the 2-d 
jet eigenvalue relationship for the wave number k is: 


tanfft (5) = 


(Ml ) 2 

f 2 ) 

\Mi) 

l/v 


(C.4) 


Equation (C.4) is used to compute the critical wave number kc and hence the wavelength This critical 
wavelength, A c is the location of the first minimum value of the interface and is more usually expressed as a 
critical location The fundamental period, T c , is found using an area weighted velocity: 



U ave = U, + U 2 (H-l) 


(C.5) 


where H denote the upper wall. These relationships are for two-dimensional geometries. Using axi-symmetric 
versions of equations (C.2)-(C.3) (Liepmann and Roshko (1957)) the eigenvalue relationship for axi-symmetric 
flows as found by Pai (1954) may be written: 
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J,(P,k) Ko (P 2 k) J Mi "| 

Jo(P,k) K,(P 2 k) l Mi) 

in terms of Bessel functions and modified Bessel functions (see Abramowitz and Stegun (1964)). Furthermore, 
a secondary stream that is stationary or incompressible yields: 

y = 2P, = 2(M]-lfi (C.7) 

which is Prandtl's (1904) result for a two-dimensional jet effluxing into a stagnant environment. The analogous 
relationships for axi-symmetric flows, Equation (C.7), for M 2 =0 are: 


A 


(C.6) 


Jo (P, k) = 0 x c = 1 .307 (M 2 , - 1 h (C.8) 

which also agrees with Prandtl’s result as quoted in Pai (1954). 


C.2.2 Slipline perturbation model 

The inlet perturbation model above is useful for small pressure disturbances but comparison with available 
experimental information is poor (Love et al. (1959)). Pack (1950), developed an alternative linear analysis for 
free jets involving a perturbation about the slipline itself. Here, for later reference. Packs model is outlined. 
Since Pack’s model is intimately related to the strained coordinate model it is appropriate to give some details 
of the analysis instead of merely quoting it. Referring to Figure C.l, the analysis perturbs the flow about the 
streamline or slipline locations. Equations (C.l)- (C.3) remain the same, except, that a location subscript, ’s’ is 
used to indicate that the perturbation has been made about the slipline values. Using a perturbation expansion 
yields: 


u d(fr 

— = J + ~r-£+... 
U Is dx 


"'-ft 


(C.9) 


The slipline perturbation method describes the flow better than the inlet perturbation because it represents 
the Prandtl-Meyer expansion in the primary stream, whereas, the inlet perturbation method merely 
approximates it. However the price is the solution of a nonlinear system of equations, the 0(1) system, 
although for the free jet case, these relationships are relatively simple and may be solved explicitly. 
Perturbation methods that retain as much nonlinearity in their lowest order expansion as practical are more 
physically complete and hence powerful relative to their less complex counterparts. 


C.2.3 Base flow slipline analysis 

The base or lowest order flow is a set of isentropic expansions (or weak compressions) in the primary and 
secondary flows. Since the flows being considered are typically under expanded, it is reasonable to expect 
expansions of the primary flow and compression of the secondary flow. Modeling of a slightly overexpanded 
primary flow is also possible because for small pressure differences, the increase in entropy for the 
compression (oblique shock) is negligible. 
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Applying the static pressure matching condition and total pressure definition (Anderson (1985)) at the 
slipline interface yields: 


Mis = Z-j[(^r)'r (l + Zy-Mb-l] (C.10) 

y-i roi 2 

For the degenerate case of a supersonic jet effluxing into a stagnant atmosphere, which was also reported 
by Pack, equation (C.10) reduces to: 


Mi = 


2 , y-1 , p, r± 

jr, [(M-Mix-f), -a 


(C.ll) 


C.3 First order slipline analysis 

With the lowest order slipline quantities available from equations (C.10) and (C.ll) the governing 
equations for the higher order, 0(e), portion of the analysis may now be derived. The equation governing the 
primary stream perturbation flow is similar to equation (C.2), except for the subscript "s M , and for an axi- 
symmetric flow: 


ga d*ti _ l 

dx 2 dr 2 r dr 


The appropriate boundary conditions for the flow are: 
• at the centerline: 


(C. 12) 


<f>, (x, 0) = finite value 


(Cl 3) 


• at the slipline defined by £,(x): 


dx 


MM i 

dx dxdr 


dx 


(C.14) 


Equation (C.14) introduces a transfer of boundary condition approximation often used in slender body theory, 
(Van Dyke (1975) and Nayfeh (1973)), and if the displacement, is small, the approximation is reasonable. It 
is this boundary condition that causes limitations as the stream pressure ratio increases. To complete the set of 
equations, the initial conditions are: 


<fi l (0,r) = 0 


d0,(O,r) 

dx 


(C. 15) 
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Since our goal is to find the slipline geometry ^(x), the use of the kinematic condition (Drazin and Reid 
(1981)) is appropriate: 


dr dx 


m=o 


(C.16) 


Equations (C.12)-(C.16) describe the primary stream flow field and the slipline location. This contrasts 
with the inlet perturbation method in which governing equations were required for both streams. The 
difference lies in that for the slipline perturbation method information about the secondary stream is contained 
in the first order solution and is passed to the higher order analysis through the coefficient p ls in equation 
(C.14), whereas, the inlet perturbation method requires explicit information about the secondary stream in the 
higher order 0(e), equation system. In addition, initial condition information is built into the slipline 
perturbation method unlike the "normal mode" method of the inlet perturbation model. 

Equations (C.12)-(C.16) have been solved using eigenfunction methods (Pack (1950)) to yield the 
potential and the displacements. Since the goal is to later obtain improvements to Pack's solution for 
unconfined flows it is necessary describe the velocity potential and displacement potentials: 


<(>(x,r) = -2p h Y J 


1 x 

^ }J , “ :Jo(2.„r)sm(X n —) 

n-1 J 1 ( A,n) r' Is 


{CM) 


and the displacement is: 


m - 4if u 



sin 2 ( 


2P ts 


J(X n ) = 0 


(C. 1 8) 


From Pack (1950), the critical location (first minimum) is at: 



(C.19) 


for axi-symmetric flows. The analogous relationship for two-dimensional flows is: 


J = (C.20) 

Equations (C.19) and (C.20) contain the pressure ratio as pj s , consequently it is expected that they are a 
improvement over the inlet perturbation solutions. The next section shows, that in spite of the improvements 
due to the slipline method, equations (C.19) and (C.20) are still limited to relatively small pressure ratios 
between streams. Additionally, equation (C.17) admits the trivial solution of: <|>](x,r) = £(x)— 0 for e=0 which is 
not experimentally or analytically valid as shown by Love et ah (1959). Love, et al. indicate that there is a 
particular solution valid for both two-dimensional and axi-symmetric flows for p 2 /pi = l which states x < /2=pj s . 
This particular solution should be recovered by equations (C.20) and (C.19) but is not. 
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C. 4 Strained coordinate model 


To overcome the limitations of Pack’s streamline perturbation analysis, namely, small pressure ratio 
restrictions and singularity for e=0, a strained coordinate slipline perturbation method has been formulated. 
The origin and full description of the strained coordinate method may be found in Van Dyke (1975), Nayfeh 
(1973) and Kevorkian and Cole (1981). A related approach is the Rayleigh-Schrodinger method described by 
Nayfeh (1973). Motivation for the use of the strained coordinates comes from the recognition, that Pack's 
analysis fails for large pressure ratios between streams that are associated with large displacement of the 
interface. This causes the boundary condition transfer analysis in equation (14) to become a weak 
approximation. Defining new coordinates that better approximate the location of the slipline velocity boundary 
condition might be expected to give a better solution. This idea is central to Lighthill's method of strained 
coordinates described by Van Dyke (1975). 

The strained coordinate method starts by considering the solution of ordinary differential equations 
resulting from a separation of variable solution of equation (C.12) with <|)(x,r)=f(x)g(r). Without the transfer of 
boundary conditions approximation of equation (C.14) the streamwise varying differential equation for f(x), is 
obtained as: 


d 2 f 

dx 2 



1 _ 

(i+mf 


f = 0 


J O (Xn) = 0 


(C.21) 


Following the slipline perturbation solution method the slipline displacement, £(x) of equation (C.18), is 
approximated as: 


#M«f.sin , rTT L ; 


2 P 


Is 





(C.22) 


Since Pack's solution is known to be singular for zero displacement, £~e=0, a stretching transformation on the 
displacement: 


4(x) 




Si 


Z Pis 


(C.23) 


is needed, where the displacement magnitude as the appropriate perturbation parameter is selected. 
Although appropriate for £ m «l, this transformation amplifies the problem of convergence of the series, see 
equation (C.9), for large displacement magnitudes. This problem of convergence, though, is precisely the type 
of situation where the method of strained coordinates is useful. Introducing new (strained) coordinates: 


X = S+X 2 (s)&+.. (C.24) 


where "s" is our new coordinate and x 2 (s) is to be determined. Similarly the basic solution is perturbed as: 

i_ 

f = fo + L 2 f,+- (C.25) 
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Equation (C.24) is used to transform derivatives such as: 


and 


d .'d 

— = 0-x 2 

dx as 


(C.26) 


£_ 

dx 2 


(l-2x\ 


' d 2 ■■ 

C2) 7 “ 

w ds 2 1 


hm ds 


(C.27) 


The expansions defined by equations (C.24) and (C.25), and the transformed derivatives are substituted into 
equation (C.21) to give to 0(1): 


d*f 9 

ds 2 


+[Y-ff 0 = o 

Pis 


(C.28) 


Solution of equation (C.28) yields the previously derived harmonic solution, which forms the basis of equation 
(Cl 7). To 0(^ m 1/2 ), equation (C.28) yields: 


d 2 f 1 , -,2 r " df o ^ * d f 0 , o . 2 /An S \ r An l 2 f 

17* l Tj frXl H ^ 2s “ ( W?TS u 


(C.29) 


To prevent the existence of solutions, which grow unphysically large, x 2 (s), is chosen so that it 
approximately eliminates the right hand side of equation (C.29). This can be considered as eliminating a 
forcing term in a harmonic system. One possible choice to remove the right hand side is to require: 

so, ha, «*» 

For large the last term in equation (C.30) may be neglected, giving x 2 (s)=s/2, then x" 2 -0 and from equation 
(C.24), and the strained coordinate "s M is given by: 


S = 


1 K L 


(C.31) 


The 0(1) solution of equation (C.21) is the same as the unstrained problem (except x=s), so that following 
Lighthill's method (as described by Van Dyke (1975)), equation (C.31) may be introduced into the solutions for 
the unstrained solutions: equations (C. 1 7) and (C.18). Further, since equation (C.31) involves straining by a 


207 



constant, the minimization of equation (CAS) remains valid to within a constant, hence equations (C.19) and 
(C.20) may be written: 

•axi-symmetric: 


y = 1-22 fi# = 1.22(1 + ^Zi)P ls (C.32) 

•two-dimensional: 

y = 2 Peff = 2(1 + l ~$j)P ls (C.33) 

Equations (C.32) and (C.33) represent the strained coordinate modifications of Pack’s slipline analysis. 
Since the problem is valid for small pressure ratios due to the stretching transformation of equation (C.32), it is 
permissible to modify the strained coordinate solutions using the particular solution referenced by Love et al. 
(1959), valid for pressure ratios precisely equal to one. Thus: 

•axi-symmetric: 


f - 1.22P*-.22P„ « 1.22(1 + L 2 (l)P h -22P„ (C.34) 

•two-dimensional: 

y = 2 Peff- Pis = 2(1 + J)P,s - P, s (C.35) 

where (3| s =(M ls 2 -l) 1 ' and ^=P i s 2 [ 1 -U i/I_T ts ]. Equations (C.34) and (C.35) are the final modified solutions for 
the critical slipline location. These relationships comprise our solution for the location of the first minimum of 
the slipline. 


C. 5 Results and comparison of models 

In this section, we compare the different models using a degenerate free jet case. The critical location and 
period of the interface, Xc and T c , are the parameters of primaiy interest in this analysis and are obtained from 
the model solutions listed in Table C.l. 
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Model Type 

Axi-symmetric 

Two-dimensional 

Inlet Perturbation (IP) 

equation (C.8) 

equation (C.7) 

Slipline Perturbation (SP) 

equation (C.l 9) 

equation (C.20) 

Strained Coordinate (SC) 

equation (C.34) 

equation (C.35) 


TABLE C.l A summary of inviscid slipline solution types. 


C.5.7 Inlet perturbation model 

Results for an axi-symmetric free jet, equation (C.8), are not presented since the model is limited to 
matched static pressure flows. Although equation (C.8) does provide an adequate comparison for sonic 
primary jet problems, it becomes poor as the Mach number increases. This poor comparison is attributed to 
both the pressure ratio between streams playing no role in the analysis and the linearization of the primary 
stream expansion. These limitations are removed using the slipline perturbation method. 


C.5.2 Slipline perturbation 

Comparison of axi-symmetric free jet results from the slipline perturbation solutions (SP) of Pack (1950) 
with the data of Love et al. (1959), are shown in Figure C.2 for several Mach numbers. Figure C.2 shows a 
good comparison with the data of Love et al. for small pressure ratios and low primary Mach numbers. In 
general, the axi-symmetric solutions for a free jet show considerable improvement over the inlet perturbation 
solutions. The improvement is attributed to the slipline perturbation containing the pressure ratio as an explicit 
parameter. However, the solution yields poor results for large pressure ratios and the error is also large for 
equal pressure ratios. These results prompted the development of the strained coordinate streamline 
perturbation model. 


C.5.3 Strained coordinates 

The strained coordinate modification of Pack's solution (SC) is compared with the axi-symmetric, free-jet 
experiment of Love et al. (1950), for several primary Mach numbers, in Figure C.2. Love et al. provides a 
semi-empirical curve fit relationship for the critical location that takes the form: 
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= 1.55[ Mi 2 0-55 p, ^<2 

2 P 2 P 2 


x± = 1.52(^) 0437 + 1.55([2 m]- 1 f 2 - 1) - 0-55 fi, + 0.5[-^—[(^-- 2)p,ji- 1] (C.36) 
2 p 2 1.55 p 2 



2 


P 2 


Equation (C.36) provides a good fit to the data over the range compared here. It should also be noted, that 
although equation (C.36) provides a useful empirical relationship it may be difficult to directly extend the 
solution to other flow problems, for example confined flows, since it does not have a completely rational basis. 



FIGURE C.2 A comparison of the strained coordinate solution with Pack’s (1950) solution and the 
experimental correlation of Love et al. (1950). 
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Figure C.2 shows that the comparison between the strained coordinate model and the experimental data is 
excellent. Inspection of the figures shows that the strained coordinate theory provides a considerable 
improvement over the simple slipline perturbation (SP) method for all primary Mach numbers and for a large 
range of pressure ratios. However, the agreement is poorer for sonic primary Mach numbers, but this is to be 
expected because the strained coordinate analysis was developed to handle large, 0(1), slipline perturbations. 
Since the slipline amplitude scales according to: 


L = Pis £ = [Mis-lfi£ 


(C. 37) 


The amplitude is expected to be relatively small for Mach numbers close to one and hence the strained 
coordinate modification is not formally appropriate. In spite of these formal limitations the strained coordinate 
model has shown good comparison (maximum of 10% relative error) for Mach numbers ranging from 
approximately sonic, Mach 1.2, to Mach 4.0. A maximum relative and root mean squared errors may be 
defined: 


err rms = 


yv — j -(— ) f 

/ j L \ 2 ^ ex per. \ ^ ' anal. J 
n v 

Y(—f 

/ j \ 'y / exper. 
i Z 


xl00% 


(C.38) 


Applying this definition to the SC method and the experimental correlation of Love, equation (C.36), yields the 
results presented in Table C.2 for a pressure ratio range of Pi/p 2 = 1.0-20.0. 


Primary Mach number 

RMS error equation (C.38) 

Maximum relative error 

1.2 

6.34 % 

9.48 % 

2.0 

1.15% 

5.44% 

3.0 

4.12% 

5.19% 

4.0 

9.03 % 

10.10% 


TABLE C.2 RMS error and maximum relative error between strained coordinate model and measurements 

of Love et al. (1959). 


Finally, it must be noted that for precisely sonic flow with pressure matching, p 2 /pi=l, the error is zero, 
but increases for larger pressure ratios as discussed previously. Therefore, the strained coordinate solution 
provides a greatly improved prediction of free jet critical locations over a full range of Mach numbers and 
pressure ratios. This concludes our discussion of the supersonic/subsonic interface problem. The use of 
these relationships in the main ejector code is described in Section 2. 


211 

















APPENDIX D 


ANALYTICAL SOLUTION OF MIXING DIFFERENTIAL EQUATIONS USING COSINE 
TRANSFORM AND METHOD OF IMAGES 


As developed in Section 2.4, the governing equations may be written for the conservation quantities in 
terms of the general linear parabolic equation (the canonical form): 


d<t> ,,d 2 <t> 

— = a(x) — - 
dx dy 2 


(D.l) 


with the boundary conditions: 


with the initial condition: 


and <[) is defined by: 


d<Mx ,0) d<f>(x,l) 

dy dy 


My) = 




0<y<hp 
h s <y< 1) 


<Mx,y) = 


( pu 2 + p' 

puH 

P“ . 


(D.2) 


(D.3) 


(D.4) 


Recall, that the turbulence model admitted the local result, that the function a(x) may be approximated 
by the linear function a(x)=a* x. The choice of the locally linear functional form is justified by the 
development of the near field turbulence model, which did indeed show a locally linear function, equation 
(2.89), e.g.: 


kxe G = 



(D.5) 


With this local approximation, it is possible to write equation (D.l) in the form: 


d<j> . d 2 
— = a x — 
dx 


dy 


(D.6) 
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wall shroud 



centerline 


FIGURE D.l Discontinuous initial conditions for the mixing problem. 


As indicated previously and repeated here in Figure D.l, the discontinuous nature of the initial condition 
(which is not even a continuous, C°, function (Marsden and Tromba (1976)), causes considerable difficulties 
both for numerical and analytical solutions. What is required is an analysis capable of dealing with extremely 
rapidly changing functions or more properly distribution. Fortunately, a solution method from classical 
analysis is available which is based upon distribution theory rather than requiring continuous functions. This is 
the method of Green's function expansions (Haberman (1983) or Weinberger (1965)). It is actually more 
useful to derive the solution using a cosine transform and the method of images, (Haberman (1983) or 
Weinberger (1965)). The ultimately desired result is: 


er f 


\ 

< (2 a) ~i x j 


-erf 


\ 

y ^hr2ri 
K (2 a 'fixj 




20 


(D.7) 


The strategy used to obtain this result is: 

• Derive a semi-infinite solution for our problem using cosine transforms (a special, semi-infinite, even 
function case of Fourier transforms). 

• Specialize this semi-infinite result for our finite problem using the method of images. 

Starting, by considering a more canonical form of the governing equation (D.6), it is necessary to 
introduce the new variable, t=l/2x 2 . To solve the semi-infinite problem then, it is necessary to consider (notice, 
also, the introduction of the dependent variable “u” which is used to signify the generality of the canonical 
form): 


du . d 2 u 
— = a — - 
dt dy 


(D.8) 


with the boundary conditions: 


’ — = 0 </)(t, y -» <*>) = 0 

dy 


(D.9) 


Initial conditions are generalized to: 
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1 0<y<hs 

o hs<y<l. 


which is shown diagramtically in Figure D.2: 



FIGURE D.2 Diagram of the canonical, step function initial conditions. 


Formally transforming equation (D.8) and solving the resulting ordinary differential equation yields: 


U(ca,t) = C(co)e~ 


Then summing (integrating) over all wave numbers, one obtains: 


w(j M)= C(a>) cos(ay)e { ° lat do) 


where the unknown function C(co) is compute using the cosine transform pair for t=0: 


C(6)) = — J^w(0,/)cos(6y)c/v 


(D.13) 


Applying our initial condition and computing the elementary integral yields the formal solution: 


o 1 

u(y,t) = — f° — s\n(a>h s )cos(coy)e~ ,ola 1 dco 
n * a) 


(D.14) 


Unfortunately, this formal solution is in terms of an infinite integral, which makes its practical evaluation 
tedious. 
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An alternative method to this situation is to return to equation (D.ll) and apply a cosine transform 
convolution theorem (Haberman (1983)): 


C- 1 (F(ffl)G(ffl)} = - )g tslfiy -y) + f O' + y)yp 
71 0 


(D.15) 


where for our problem F(a>) is the exponential term in equation (D.ll) and G(ca) is the step function initial 
condition. Fortunately, the inversion of the exponential term is tabulated: 


C" 1 



e 


-o) 2 /Aa 


= 6“” 


(D.16) 


and the step function merely carries through to yield: 


u(y,t) = C-'{U{co,t)} 


1 




Tta t 



0->T | 

4 at J 



(D.17) 


The step function has been used to evaluate the infinite integral. From this point, the integrals in equation 
(D.17) are evaluated using straightforward substitution, which yields the expected error function (Abramowitz 
and Stegun (1964)) form: 
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f \ 


f \ 


1 

erf 

(y + K) 
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-erf 

iy-K) 
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U CM) = - 


2 







Equation (D.18) has been derived to satisfy the unit step function initial condition. By elementary 
substitution one may modify this condition to apply for the physical initial condition required by equation 
(D.3): 




f \ 


r \ 


u sem i- inf (y,t) = ^(^ 10 -^ 20 ) 

erf 

(y+h.) 
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-erf 

(y - K) 
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(4a‘t)> J 


l( 4a ‘t)L 



(D.19) 


Equation (D.19) satisfies all the conditions (initial, governing equation and zero boundary condition), but 
will fail to respect the wall Nuemann condition , du/5y(y=l)=0. To satisfy this condition, it is necessary to 
invoke the method of images. Literally, it is necessary to place a family of image sources such that the new 
condition is respected. Equation (D.19) may be thought of as two sources of equal strength placed on either 
side of the y axis origin. It is possible to extend this idea (basically a symmetry operation, which uses linear 
superposition) to satisfy the entire finite condition. Following figure D.3: 
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images 


physical 


images 


y=-2 y=-l y=0 y=l y=2 

FIGURE D.3 Image “sources” used to satisfy the y=l, Nuemann condition. 

and changing variables, t=l/2x 2 , yields: 
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(D.20) 


Which is the expected solution, equation (D.7). 

Now, since the project uses this solution in a combined analytical/numerical solution scheme, it is useful to 
include the associated partial derivatives: 
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and 
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(D.22) 


Either by direct computation or referring to the governing equation, equation (D.6) it is clear that: 
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(D.23) 
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Which verifies that this solution does indeed satisfy the governing equation. This concludes the discussion of 
the method used to compute the analytical portion of our equation system. 
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APPENDIX E 


AN APPROXIMATE DERIVATION OF THE “VIGNERON” PARAMETER, co 


In this appendix a value for the Vigneron parameter, co, (Vigneron et al. (1978)) is derived by analyzing an 
approximate family of equations. Though approximate, this development provides a direct connection to the 
governing equations used in this project. More compete derivations are possible (see Anderson et al. (1984)), 
but the formulation described in this appendix, exploits the basic form of our parabolic equations. Consider the 
streamwise momentum equation: 


^-(Pu+cop) 

ox 


dyKU dy ) 


(E.l) 


Now, further assuming that the mass and energy fluxes are both constant: pu=constant and puH=constant (this 
assumption could be relaxed in favor of the full transport equations for energy and mass, but the eigenvalue 
problem that would be need to be resolved is much more complex) and rewriting state: 


p = pRT 


(E.2) 


in terms of energy yields: 


P = 


Lzl 

r 


pH-^pu 2 


(E.3) 


At this point, it is desirable to introduce the co parameter as a “marker” on the static pressure, since it is the 
static pressure modeling which will control the eigenvalue consistency. Additionally, it is possible to factor the 
constant mass flux to yield: 
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dx 



oKr - 1) 
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co(y - 1 )H 


uj 
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V U 


Re-writing, equation (E.4) in terms of u only, one obtains: 


(i _ *>(r-P )xj _ gj 

y y u 2 ) dx 



(E.4) 


(E.5) 


For a well posed diffusion problem (i.e. one with a positive eigenvalues or alternatively a positive 
“viscosity”) the term in the brackets of equation (E.5) must be greater than zero: 



o>(r - i) j M _ <o(y- 1) h ' | 
y y u 2 ) 


>0 


(E.6) 
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The similarity of this analysis with the model equation discussed in Section 2.4.3.3 is, of course, not 
coincidence. Inequality (E.6) may be simplified by recognizing that: 


7/ 1 1 

u 2 ~ (y - 1 )M 2 2 


(E.7) 


Combining equations (E.6) and (E.7) and solving the inequality for the threshold value of ca that provides a 
physical eigenvalue, yields: 


yM 2 

CO = r 

1 + (/ - 1 )M 2 


(E.8) 


Which is, of course, the expected value for the Vigneron parameter. 
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